text:
function i4_bset_test ( )

%*****************************************************************************80
%
%% I4_BSET_TEST tests I4_BSET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  test_num = 2;
  i4_test = [ 101, -31 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03225\n' );
  fprintf ( 1,'  I4_BSET sets a given bit to 0.\n' );

  for test = 1 : test_num

    i4 = i4_test(test);

    ivec = i4_to_bvec ( i4, 32 );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Working on I4 = %d\n', i4 );
    fprintf ( 1, '\n' );
    fprintf ( 1, '       Pos     Digit       I4_BSET\n' );
    fprintf ( 1, '\n' );

    for pos = 0 : 31

      j1 = i4_bset ( i4, pos );

      fprintf ( 1, '  %8d  %8d  %8d\n', pos, ivec(pos+1), j1 );

    end

  end

  return
end


text:
function value = r4_aint ( x )

%****************************************************************************80
%
%% R4_AINT truncates an R4 argument to an integer.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2011
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the truncated version of X.
%
  if ( x < 0.0 )
    value = - floor ( abs ( x ) );
  else
    value =   floor ( abs ( x ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function imgOut = pyrGaussGenAux(img)
%
%
%        imgOut = pyrGaussGenAux(img)
%
%
%        Input:
%           -img: an image
%
%        Output:
%           -imgOut: filtered img using 5x5 Gaussian filter
%           -imgB: downsampled img at half the size
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

%5x5 Gaussian Kernel
kernel = [1, 4, 6, 4, 1];
mtx = kernel' * kernel;
mtx = mtx / sum(mtx(:));

%Convolution
imgB = imfilter(img, mtx, 'replicate');

%Downsampling
[r, c] = size(img);
imgOut = imgB(1:2:r, 1:2:c); %imresize(imgB, 0.5, 'bilinear');

end

author: banterle


text:
function [ y ] = back_trans( x,waveS,waveletFilterName )
%FFT3D Summary of this function goes here
%   Detailed explanation goes here

if isreal(x{1,1})
    y = cell(1,size(x,2));
    for j = 1:size(x,2)
        y{1,j} = waverec2(x{1,j},waveS,waveletFilterName);
    end;
else
    y = cell(1,2*size(x,2));
    for j = 1:size(x,2)
        y{1,j} = waverec2(real(x{1,j}),waveS,waveletFilterName);
        y{1,j+size(x,2)} = waverec2(imag(x{1,j}),waveS,waveletFilterName);
    end;
end;
end

author: thomaskuestner


text:
function linplus_test43 ( )

%*****************************************************************************80
%
%% TEST43 tests R8LT_MXV, R8LT_SL, R8LT_VXM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST43\n' );
  fprintf ( 1, '  For a matrix in lower triangular storage,\n' );
  fprintf ( 1, '  R8LT_SL solves systems;\n' );
  fprintf ( 1, '  R8LT_MXV computes matrix-vector products;\n' );
  fprintf ( 1, '  R8LT_VXM computes vector-matrix products;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );

  for i = 1 : n
    for j = 1 : n
      if ( j <= i )
        a(i,j) = j;
      else
        a(i,j) = 0.0;
      end
    end
  end

  r8lt_print ( n, n, a, '  The lower triangular matrix:' );

  for job = 0 : 1
%
%  Set the desired solution.
%
    x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
    if ( job == 0 )
      b = r8lt_mxv ( n, n, a, x );
    else
      b = r8lt_vxm ( n, n, a, x );
    end
%
%  Solve the linear system.
%
    x = r8lt_sl ( n, a, b, job );

    if ( job == 0 )
      r8vec_print ( n, x, '  Solution:' );
    else
      r8vec_print ( n, x, '  Solution to the transposed system:' );
    end

  end

  return
end


text:
function results = vl_test_fisher(varargin)
% VL_TEST_FISHER
vl_test_init ;

function s =  setup()
randn('state',0) ;
dimension = 5 ;
numData = 21 ;
numComponents = 3 ;
s.x = randn(dimension,numData) ;
s.mu = randn(dimension,numComponents) ;
s.sigma2 = ones(dimension,numComponents) ;
s.prior = ones(1,numComponents) ;
s.prior = s.prior / sum(s.prior) ;

function test_basic(s)
phi_ = simple_fisher(s.x, s.mu, s.sigma2, s.prior) ;
phi = vl_fisher(s.x, s.mu, s.sigma2, s.prior) ;
vl_assert_almost_equal(phi, phi_, 1e-10) ;

function test_norm(s)
phi_ = simple_fisher(s.x, s.mu, s.sigma2, s.prior) ;
phi_ = phi_ / norm(phi_) ;
phi = vl_fisher(s.x, s.mu, s.sigma2, s.prior, 'normalized') ;
vl_assert_almost_equal(phi, phi_, 1e-10) ;

function test_sqrt(s)
phi_ = simple_fisher(s.x, s.mu, s.sigma2, s.prior) ;
phi_ = sign(phi_) .* sqrt(abs(phi_)) ;
phi = vl_fisher(s.x, s.mu, s.sigma2, s.prior, 'squareroot') ;
vl_assert_almost_equal(phi, phi_, 1e-10) ;

function test_improved(s)
phi_ = simple_fisher(s.x, s.mu, s.sigma2, s.prior) ;
phi_ = sign(phi_) .* sqrt(abs(phi_)) ;
phi_ = phi_ / norm(phi_) ;
phi = vl_fisher(s.x, s.mu, s.sigma2, s.prior, 'improved') ;
vl_assert_almost_equal(phi, phi_, 1e-10) ;

function enc = simple_fisher(x, mu, sigma2, pri)
sigma = sqrt(sigma2) ;
for i = 1:size(mu,2)
  delta{i} = bsxfun(@times, bsxfun(@minus, x, mu(:,i)), 1./sigma(:,i)) ;
  q(i,:) = log(pri(i)) - 0.5 * log(sigma2(i)) - 0.5 * sum(delta{i}.^2,1) ;
end
q = exp(bsxfun(@minus, q, max(q,[],1))) ;
q = bsxfun(@times, q, 1 ./ sum(q,1)) ;
n = size(x,2) ;
for i = 1:size(mu,2)
  u{i} = delta{i} * q(i,:)' / n / sqrt(pri(i)) ;
  v{i} = (delta{i}.^2 - 1) * q(i,:)' / n / sqrt(2*pri(i)) ;
end
enc = cat(1, u{:}, v{:}) ;

author: jianxiongxiao

repo: ProfXkit


text:
function tf = strlexle(a, b)
%STRLEXLE Lexicographic less than or equal to.
%
%   STRLEXLE(A, B) returns 1 if A is lexicographically less than or equal to B,
%   and 0 otherwise.
%
%   This is a MATLAB version of the Perl `le' operator.
%
%   See also LE.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-10-13 11:15:58 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   tf = strlexcmp(a, b) <= 0;

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function userInterrupt(standalone)
global procStreamStyle
global testidx

UnitTests_Init(standalone, false);
testidx=[];
procStreamStyle=[];

author: BUNPC


text:
function d = mahalanobis(varargin)
%MAHALANOBIS Computes the Mahalanobis distance.
%   D = MAHALANOBIS(Y, X) computes the Mahalanobis distance between
%   each vector in Y to the mean (centroid) of the vectors in X, and
%   outputs the result in vector D, whose length is size(Y, 1).  The
%   vectors in X and Y are assumed to be organized as rows.  The
%   input data can be real of complex. The outputs are real
%   quantities.
%
%   D = MAHALANOBIS(Y, CX, MX) computes the Mahalanobis distance
%   between each vector in Y and the given mean vector, MX. The
%   results are output in vector D, whose length is size(Y, 1).  The
%   vectors in Y are assumed to be organized as the rows of this
%   array. The input data can be real or complex. The outputs are
%   real quantities. In addition to the mean vector MX, the
%   covariance matrix CX of a population of vectors X also must be
%   provided. Use function COVMATRIX (Section 11.5) to compute MX and
%   CX.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/10/26 23:19:44 $

% Reference: Acklam, P. J. [2002]. "MATLAB Array Manipulation Tips
% and Tricks." Available at
%     home.online.no/~pjacklam/matlab/doc/mtt/index.html
% or at
%     www.prenhall.com/gonzalezwoodseddins

param = varargin; % Keep in mind that param is a cell array.
Y = param{1};
ny = size(Y, 1); % Number of vectors in Y.

if length(param) == 2
   X = param{2};
   % Compute the mean vector and covariance matrix of the vectors
   % in X.
   [Cx, mx] = covmatrix(X);
elseif length(param) == 3 % Cov. matrix and mean vector provided.
   Cx = param{2};
   mx = param{3};
else
   error('Wrong number of inputs.')
end
mx = mx(:)'; % Make sure that mx is a row vector.

% Subtract the mean vector from each vector in Y.
Yc = Y - mx(ones(ny, 1), :);

% Compute the Mahalanobis distances.
d = real(sum(Yc/Cx.*conj(Yc), 2));

repo: weiyanmin


text:
function [Hjorth] = eeg_lap_hjorth(voltage,X,Y)

% eeg_lap_hjorth - 2D Laplacian of Potential at XY
%
% Useage: [Hjorth] = eeg_lap_hjorth(voltage [,X,Y])
%
% The Hjorth nearest neighbour, finite difference Laplacian.
% For a continuous approximation of the Laplacian, use the
% 'eeg_lap' function.
%
% This routine simply calls the del2 matlab command, which requires
% 'voltage' to be a rectangular matrix.  See 'help del2', esp:
%
% L = DEL2(U,HX,HY) when U is 2-D, uses the spacing specified by HX
%     and HY. If HX is a scalar, it gives the spacing between points in
%     the x-direction. If HX is a vector, it must be of length SIZE(U,2)
%     and specifies the x-coordinates of the points.  Similarly, if HY
%     is a scalar, it gives the spacing between points in the
%     y-direction. If HY is a vector, it must be of length SIZE(U,1) and
%     specifies the y-coordinates of the points.
%
% For example:
%
% [x,y] = meshgrid(-10:.5:10);
% z = (x.^2).*(y.^2); % simulate monotonic potential
% del2z = 4*del2(z);  % calculate Hjorth laplacian
% figure('name','potential vs laplacian','numbertitle','off','position',[500 10 512 512]);
% subplot(2,1,1); surf(x,y,z);
% title('potential'); shading interp; colorbar; rotate3d; axis tight
% subplot(2,1,2); surf(x,y,del2z),
% title('Hjorth laplacian'); shading interp; colorbar; rotate3d; axis tight
%
% refs:  Hjorth B (1975).  An on-line transformation of EEG scalp
%          potentials into orthogonal source derivations.
%          Electroencephalography & Clinical Neurophysiology, 39: 526-530.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:51 $

% Licence:  GNU GPL, no implied or express warranties
% History:  06/01, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('X','var') X = 1; end
if ~exist('Y','var') Y = 1; end

Hjorth = 4*del2(voltage,X,Y);

return

author: PatternRecognition

repo: OpenBMI


text:
%*****************************************************************************
% DSDP5:  Dual Scaling Algorithm for Positive Semidefinite Programming
% Copyright (c) 2002 by
% S. J. Benson, Y. Ye
% Last modified: 20 January 2004
%*****************************************************************************
%
% Converts data from DSDP4 format to DSDP5 format, calls solver, and
% and converts solution to DSDP4 format.
%
%
% > DSDP(A,C,b) attempts to solve the positive semidefinite program
%      MINIMIZE trace(C*X) SUCH THAT trace(A_i*X) = b_i, i=1,...,m and X >= 0
%      using a dual scaling algorithm.  For a problem with p blocks and m
%      constraints, A is a p x m cell array and C is a p x 1 cell array.
%      One block may contain LP variables, and the cells corresponding to
%      this block should be a one dimensional array.  All other cells
%      must contain a square, symmetric, real valued matrix.  The third
%      argument b is a dense column vector of length m.
%
% > DSDP(A,C,b,OPTIONS,y0) specifies an initial dual vector y0.
%
% > [STAT,y,X] = DSDP() returns a structure containing relevant statistics,
%                and approximate dual and primal solutions.
%
%*****************************************************************************

function [STAT,y,XX] = dsdp4(A,C,b,OPTIONS,y0);

  p=length(C);
  m=length(b);

  AC=cell(p,3);
  for j=1:p,
     [n1,n2]=size(C{j});
     if (n1==1 | n2==1)
       AAC=sparse(n1*n2,m+1);
       for i=1:m, AAC=[AAC sparse(A{j,i})']; end;
       AAC=[AAC sparse(C{j})'];
       AC{j,1}='LP';
       AC{j,2}=length(C{j});
       AC{j,3}=AAC;
     else
       AAC=sparse(n1*(n1+1)/2,m+1);
       for i=1:m, AAC=[AAC dvec(A{j,i})]; end;
       AAC=[AAC dvec(C{j})];
       AC{j,1}='SDP';
       AC{j,2}=size(C{j},1);
       AC{j,3}=AAC;
     end;
  end;

  [STAT,y,X]=dsdp(b,AC,OPTIONS,y0);

  XX=cell(p,1);
  for j=1:p,
     [n1,n2]=size(C{j})
     if (n1==1 | n2==1)
       XX{j}=X{j};
     else
       XX{j}=dmat(X{j});
     end;
  end;

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [ compSymbolList, compNameList  ] = getDefaultCompartments( )
%GETDEFAULTCOMPARTMENTS returns the default compartment Symbols and default Compartments
%
%   USAGE:  [ compSymbolList, compNameList  ] = getDefaultCompartments( )
%
% OUTPUT:
%    compSymbolList:       Default symbols of compartments
%    compNameList:         Names of the default compartments.
%
% .. Authors: Thomas Pfau May 2017
    compSymbolList = {'c','m','v','x','e','t','g','r','n','p','l','u'};
    compNameList = {'Cytoplasm','Mitochondrion','Vacuole','Peroxisome','Extracellular','Pool','Golgi','Endoplasmic_reticulum','Nucleus','Periplasm','Lysosome','Unknown'};

end

author: opencobra

repo: cobratoolbox


text:
function F = times(F, G)
%.*   Pointwise multiplication of two CHEBFUN3V objects.
%   F.*G computes pointwise, componentwise multiplication of two
%   CHEBFUN3V objects, a CHEBFUN3V and a CHEBFUN3, or a CHEBFUN3V
%   and a double.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty(F) || isempty(G) )
    F = chebfun3v();
    return
end

if ( ~isa(F, 'chebfun3v') )
    F = times(G, F);
    return
end

nF = F.nComponents;

if ( isa(G, 'double') )             % CHEBFUN3V .* double
    if ( numel(G) == 1 )            % CHEBFUN3V .* scalar
        scalar = G;
        for jj = 1:nF
            F.components{jj} = times(F.components{jj}, scalar);
        end
    elseif ( (size(G, 1) == nF) || ( F.isTransposed && (size(G, 2) == nF) ) )
                                    % CHEBFUN3V .* vector
        for jj = 1 : nF
            F.components{jj} = times(F.components{jj}, G(jj));
        end
    else
        error('CHEBFUN:CHEBFUN3V:times:double', ...
            'CHEBFUN3V and double size mismatch.');
    end

elseif ( isa(G, 'chebfun3v') )      % CHEBFUN3V .* CHEBFUN3V
    nG = G.nComponents;
    if ( nF ~= nG )
         error('CHEBFUN:CHEBFUN3V:times:numel', ...
             'CHEBFUN3V components mismatch.');
    end
    for jj = 1:nF
        F.components{jj} = times(F.components{jj}, G.components{jj});
    end

elseif ( isa(G, 'chebfun3') )       % CHEBFUN3V .* CHEBFUN3
    for jj = 1 : nF
        F.components{jj} = times(F.components{jj}, G);
    end

else  % error
    error('CHEBFUN:CHEBFUN3V:times:inputs', 'Unrecognized input arguments.');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function a = gt( x, y )

%Disciplined convex programming information for GT (>):
%   The right-hand side of a less-than constraint must be convex. The
%   left-hand side must be concave. Of course, real constant and affine
%   expressions are both convex and concave and can be used on either
%   side as well.
%
%Disciplined geometric programming information for GT (>):
%   The right-hand side of a less-than constraint must be log-convex---
%   including positive constants, monomials, posynomials, generalized
%   posynomials, and products thereof. The left-hand side must be log-
%   concave---including positive constants, monomials, reciprocals of
%   log-convex expressions, and products thereof.
%
%Note that CVX does not distinguish between strict greater-than (>) and
%greater-than-or-equal (<=) constraints; they are treated identically.
%Feasible interior-point solvers tend to return points which satisfy
%strict inequality, but not all solvers do.

b = newcnstr( evalin( 'caller', 'cvx_problem', '[]' ), x, y, '>' );
if nargout, a = b; end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [M, RHS, Mx, My, Mz, RHSx, RHSy, RHSz] = ...
    convectionTvdTermCylindrical3D(u, phi, FL)
% This function uses the TVD scheme to discretize a 3D
% convection term in the form \grad (u \phi) where u is a face vactor
% It also returns the x and y parts of the matrix of coefficient.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% extract data from the mesh structure
Nr = u.domain.dims(1);
Ntetta = u.domain.dims(2);
Nz = u.domain.dims(3);
G=reshape((1:(Nr+2)*(Ntetta+2)*(Nz+2)), Nr+2, Ntetta+2, Nz+2);
DRp = repmat(u.domain.cellsize.x(2:end-1), 1, Ntetta, Nz);
DTHETAp = repmat(u.domain.cellsize.y(2:end-1)', Nr, 1, Nz);
DZ = zeros(1,1,Nz+2);
DZ(1,1,:) = u.domain.cellsize.z;
DZp=repmat(DZ(1,1,2:end-1), Nr, Ntetta, 1);
rp = repmat(u.domain.cellcenters.x, 1, Ntetta, Nz);
rf = repmat(u.domain.facecenters.x, 1, Ntetta, Nz);
dx=repmat(0.5*(u.domain.cellsize.x(1:end-1)+u.domain.cellsize.x(2:end)), 1, Ntetta, Nz);
dy=repmat(0.5*(u.domain.cellsize.y(1:end-1)+u.domain.cellsize.y(2:end))', Nr, 1, Nz);
dz=zeros(1, 1, Nz+1);
dz(1,1,:)=0.5*(u.domain.cellsize.z(1:end-1)+u.domain.cellsize.z(2:end));
dz=repmat(dz, Nr, Ntetta, 1);
psiX_p = zeros(Nr+1,Ntetta,Nz);
psiX_m = zeros(Nr+1,Ntetta,Nz);
psiY_p = zeros(Nr,Ntetta+1,Nz);
psiY_m = zeros(Nr,Ntetta+1,Nz);
psiZ_p = zeros(Nr,Ntetta,Nz+1);
psiZ_m = zeros(Nr,Ntetta,Nz+1);

% define the vectors to stores the sparse matrix data
iix = zeros(3*(Nr+2)*(Ntetta+2)*(Nz+2),1);
jjx = zeros(3*(Nr+2)*(Ntetta+2)*(Nz+2),1);
sx = zeros(3*(Nr+2)*(Ntetta+2)*(Nz+2),1);
iiy = zeros(3*(Nr+2)*(Ntetta+2)*(Nz+2),1);
jjy = zeros(3*(Nr+2)*(Ntetta+2)*(Nz+2),1);
sy = zeros(3*(Nr+2)*(Ntetta+2)*(Nz+2),1);
iiz = zeros(3*(Nr+2)*(Ntetta+2)*(Nz+2),1);
jjz = zeros(3*(Nr+2)*(Ntetta+2)*(Nz+2),1);
sz = zeros(3*(Nr+2)*(Ntetta+2)*(Nz+2),1);
mnx = Nr*Ntetta*Nz;	mny = Nr*Ntetta*Nz;   mnz = Nr*Ntetta*Nz;

% extract the velocity data
% note: size(ux) = [1:m+1, 1:n] and size(uy) = [1:m, 1:n+1]
ux = u.xvalue;
uy = u.yvalue;
uz = u.zvalue;

% calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
% x direction
dphiX_p = (phi.value(2:Nr+2, 2:Ntetta+1, 2:Nz+1)-phi.value(1:Nr+1, 2:Ntetta+1, 2:Nz+1))./dx;
rX_p = dphiX_p(1:end-1,:,:)./fsign(dphiX_p(2:end,:,:));
psiX_p(2:Nr+1,:,:) = 0.5*FL(rX_p).* ...
    (phi.value(3:Nr+2,2:Ntetta+1,2:Nz+1)-phi.value(2:Nr+1,2:Ntetta+1,2:Nz+1));
psiX_p(1,:,:) = 0; % left boundary
% y direction
dphiY_p = (phi.value(2:Nr+1, 2:Ntetta+2, 2:Nz+1)-phi.value(2:Nr+1, 1:Ntetta+1, 2:Nz+1))./dy;
rY_p = dphiY_p(:,1:end-1,:)./fsign(dphiY_p(:,2:end,:));
psiY_p(:,2:Ntetta+1,:) = 0.5*FL(rY_p).* ...
    (phi.value(2:Nr+1,3:Ntetta+2,2:Nz+1)-phi.value(2:Nr+1, 2:Ntetta+1,2:Nz+1));
psiY_p(:,1,:) = 0; % Bottom boundary
% z direction
dphiZ_p = (phi.value(2:Nr+1, 2:Ntetta+1, 2:Nz+2)-phi.value(2:Nr+1, 2:Ntetta+1, 1:Nz+1))./dz;
rZ_p = dphiZ_p(:,:,1:end-1)./fsign(dphiZ_p(:,:,2:end));
psiZ_p(:,:,2:Nz+1) = 0.5*FL(rZ_p).* ...
    (phi.value(2:Nr+1,2:Ntetta+1,3:Nz+2)-phi.value(2:Nr+1,2:Ntetta+1,2:Nz+1));
psiZ_p(:,:,1) = 0; % Back boundary

% calculate the upstream to downstream gradient ratios for u<0 (- ratio)
% x direction
rX_m = dphiX_p(2:end,:,:)./fsign(dphiX_p(1:end-1,:,:));
psiX_m(1:Nr,:,:) = 0.5*FL(rX_m).* ...
    (phi.value(1:Nr, 2:Ntetta+1, 2:Nz+1)-phi.value(2:Nr+1, 2:Ntetta+1, 2:Nz+1));
psiX_m(Nr+1,:,:) = 0; % right boundary
% y direction
rY_m = dphiY_p(:,2:end,:)./fsign(dphiY_p(:,1:end-1,:));
psiY_m(:,1:Ntetta,:) = 0.5*FL(rY_m).* ...
    (phi.value(2:Nr+1,1:Ntetta,2:Nz+1)-phi.value(2:Nr+1,2:Ntetta+1,2:Nz+1));
psiY_m(:,Ntetta+1,:) = 0; % top boundary
% z direction
rZ_m = dphiZ_p(:,:,2:end)./fsign(dphiZ_p(:,:,1:end-1));
psiZ_m(:,:,1:Nz) = 0.5*FL(rZ_m).* ...
    (phi.value(2:Nr+1,2:Ntetta+1,1:Nz)-phi.value(2:Nr+1,2:Ntetta+1,2:Nz+1));
psiZ_m(:,:,Nz+1) = 0; % front boundary
% reassign the east, west, north, and south velocity vectors for the
% code readability
ue = ux(2:Nr+1,:,:);		uw = ux(1:Nr,:,:);
vn = uy(:,2:Ntetta+1,:);     vs = uy(:,1:Ntetta,:);
wf = uz(:,:,2:Nz+1);     wb = uz(:,:,1:Nz);
re = rf(2:Nr+1,:,:);         rw = rf(1:Nr,:,:);

% find the velocity direction for the upwind scheme
ue_min = min(ue,0);	ue_max = max(ue,0);
uw_min = min(uw,0);	uw_max = max(uw,0);
vn_min = min(vn,0);	vn_max = max(vn,0);
vs_min = min(vs,0);	vs_max = max(vs,0);
wf_min = min(wf,0);	wf_max = max(wf,0);
wb_min = min(wb,0);	wb_max = max(wb,0);

% calculate the coefficients for the internal cells
AE = re.*ue_min./(DRp.*rp);
AW = -rw.*uw_max./(DRp.*rp);
AN = vn_min./(DTHETAp.*rp);
AS = -vs_max./(DTHETAp.*rp);
AF = wf_min./DZp;
AB = -wb_max./DZp;
APx = (re.*ue_max-rw.*uw_min)./(DRp.*rp);
APy = (vn_max-vs_min)./(DTHETAp.*rp);
APz = (wf_max-wb_min)./DZp;

% Also correct for the boundary cells (not the ghost cells)
% Left boundary:
APx(1,:,:) = APx(1,:,:)-rw(1,:,:).*uw_max(1,:,:)./(2*rp(1,:,:)*DRp(1));   AW(1,:,:) = AW(1,:,:)/2;
% Right boundary:
AE(end,:,:) = AE(end,:,:)/2;    APx(end,:,:) = APx(end,:,:)+re(end,:,:).*ue_min(end,:,:)./(2*DRp(end)*rp(end,:,:));
% Bottom boundary:
APy(:,1,:) = APy(:,1,:)-vs_max(:,1,:)./(2*DTHETAp(1)*rp(:,1,:));   AS(:,1,:) = AS(:,1,:)/2;
% Top boundary:
AN(:,end,:) = AN(:,end,:)/2;    APy(:,end,:) = APy(:,end,:)+vn_min(:,end,:)./(2*DTHETAp(end)*rp(:,end,:));
% Back boundary:
APz(:,:,1) = APz(:,:,1)-wb_max(:,:,1)/(2*DZp(1));   AB(:,:,1) = AB(:,:,1)/2;
% Front boundary:
AF(:,:,end) = AF(:,:,end)/2;    APz(:,:,end) = APz(:,:,end) + wf_min(:,:,end)/(2*DZp(end));

AE = reshape(AE,mnx,1);
AW = reshape(AW,mnx,1);
AN = reshape(AN,mny,1);
AS = reshape(AS,mny,1);
AF = reshape(AF,mnz,1);
AB = reshape(AB,mnz,1);
APx = reshape(APx,mnx,1);
APy = reshape(APy,mny,1);
APz = reshape(APz,mnz,1);

% build the sparse matrix based on the numbering system
rowx_index = reshape(G(2:Nr+1,2:Ntetta+1,2:Nz+1),mnx,1); % main diagonal x
iix(1:3*mnx) = repmat(rowx_index,3,1);
rowy_index = reshape(G(2:Nr+1,2:Ntetta+1,2:Nz+1),mny,1); % main diagonal y
iiy(1:3*mny) = repmat(rowy_index,3,1);
rowz_index = reshape(G(2:Nr+1,2:Ntetta+1,2:Nz+1),mnz,1); % main diagonal z
iiz(1:3*mnz) = repmat(rowz_index,3,1);
jjx(1:3*mnx) = [reshape(G(1:Nr,2:Ntetta+1,2:Nz+1),mnx,1); reshape(G(2:Nr+1,2:Ntetta+1,2:Nz+1),mnx,1); reshape(G(3:Nr+2,2:Ntetta+1,2:Nz+1),mnx,1)];
jjy(1:3*mny) = [reshape(G(2:Nr+1,1:Ntetta,2:Nz+1),mny,1); reshape(G(2:Nr+1,2:Ntetta+1,2:Nz+1),mny,1); reshape(G(2:Nr+1,3:Ntetta+2,2:Nz+1),mny,1)];
jjz(1:3*mnz) = [reshape(G(2:Nr+1,2:Ntetta+1,1:Nz),mnz,1); reshape(G(2:Nr+1,2:Ntetta+1,2:Nz+1),mnz,1); reshape(G(2:Nr+1,2:Ntetta+1,3:Nz+2),mnz,1)];
sx(1:3*mnx) = [AW; APx; AE];
sy(1:3*mny) = [AS; APy; AN];
sz(1:3*mnz) = [AB; APz; AF];

% calculate the TVD correction term
div_x = -(1./(DRp.*rp)).*(re.*(ue_max.*psiX_p(2:Nr+1,:,:)+ue_min.*psiX_m(2:Nr+1,:,:))- ...
              rw.*(uw_max.*psiX_p(1:Nr,:,:)+uw_min.*psiX_m(1:Nr,:,:)));
div_y = -(1./(DTHETAp.*rp)).*((vn_max.*psiY_p(:,2:Ntetta+1,:)+vn_min.*psiY_m(:,2:Ntetta+1,:))- ...
              (vs_max.*psiY_p(:,1:Ntetta,:)+vs_min.*psiY_m(:,1:Ntetta,:)));
div_z = -(1./DZp).*((wf_max.*psiZ_p(:,:,2:Nz+1)+wf_min.*psiZ_m(:,:,2:Nz+1))- ...
              (wb_max.*psiZ_p(:,:,1:Nz)+wb_min.*psiZ_m(:,:,1:Nz)));

% define the RHS Vector
RHS = zeros((Nr+2)*(Ntetta+2)*(Nz+2),1);
RHSx = zeros((Nr+2)*(Ntetta+2)*(Nz+2),1);
RHSy = zeros((Nr+2)*(Ntetta+2)*(Nz+2),1);
RHSz = zeros((Nr+2)*(Ntetta+2)*(Nz+2),1);

% assign the values of the RHS vector
row_index = rowx_index;
RHS(row_index) = reshape(div_x+div_y+div_z,Nr*Ntetta*Nz,1);
RHSx(rowx_index) = reshape(div_x,Nr*Ntetta*Nz,1);
RHSy(rowy_index) = reshape(div_y,Nr*Ntetta*Nz,1);
RHSz(rowz_index) = reshape(div_z,Nr*Ntetta*Nz,1);

% build the sparse matrix
kx = 3*mnx;
ky = 3*mny;
kz = 3*mnz;
Mx = sparse(iix(1:kx), jjx(1:kx), sx(1:kx), (Nr+2)*(Ntetta+2)*(Nz+2), (Nr+2)*(Ntetta+2)*(Nz+2));
My = sparse(iiy(1:ky), jjy(1:ky), sy(1:ky), (Nr+2)*(Ntetta+2)*(Nz+2), (Nr+2)*(Ntetta+2)*(Nz+2));
Mz = sparse(iiz(1:kz), jjz(1:kz), sz(1:kz), (Nr+2)*(Ntetta+2)*(Nz+2), (Nr+2)*(Ntetta+2)*(Nz+2));
M = Mx + My + Mz;

end

function phi_out = fsign(phi_in)
% This function checks the value of phi_in and assigns an eps value to the
% elements that are less than or equal to zero, while keeping the signs of
% the nonzero elements
    phi_out = (abs(phi_in)>=eps).*phi_in+eps*(phi_in==0)+eps*(abs(phi_in)<eps).*sign(phi_in);
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [d] =  testing(a,d)

  %name=[]; if ~isempty(d) name=d.name; end; % record old name
  d=generate(a);
  %d.name=[name ' -> ' d.name];

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function outsig=frsyn(F,insig);
%FRSYN  Frame synthesis operator
%   Usage: f=frsyn(F,c);
%
%   `f=frsyn(F,c)` constructs a signal *f* from the frame coefficients *c*
%   using the frame *F*. The frame object *F* must have been created using
%   |frame|.
%
%   Examples:
%   ---------
%
%   In the following example a signal *f* is constructed through the frame
%   synthesis operator using a Gabor frame. The coefficients associated with
%   this Gabor expansion are contained in an identity matrix. The identity
%   matrix corresponds to a diagonal in the time-frequency plane, that is,
%   one atom at each time position with increasing frequency.:::
%
%      a = 10;
%      M = 40;
%
%      F = frame('dgt', 'gauss', a, M);
%
%      c = framenative2coef(F, eye(40));
%
%      f = frsyn(F, c);
%
%   See also: frame, frana, plotframe

complainif_notenoughargs(nargin,2,'FRSYN');
complainif_notvalidframeobj(F,'FRSYN');

L=framelengthcoef(F,size(insig,1));

F=frameaccel(F,L);

outsig=F.frsyn(insig);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% [INPUT]
% file = A string representing the full path to the Excel spreadsheet containing the dataset.
% version = A string representing the version of the dataset.
% date_format_base = A string representing the base date format used in the Excel spreadsheet for all elements except balance sheet ones (optional, default='dd/mm/yyyy').
% date_format_balance = A string representing the date format used in the Excel spreadsheet for balance sheet elements (optional, default='QQ yyyy').
% shares_type = A string representing the type of data included in the Shares sheet (optional, default='P'):
%   - 'P' for prices;
%   - 'R' for returns.
% crises_type = A string representing the type of data included in the Crises sheet (optional, default='R'):
%   - 'E' for events;
%   - 'R' for time ranges.
% distress_threshold = A float [0.05,0.75] representing the threshold at which distress statuses are detected (optional, default=0.05).
%
% [OUTPUT]
% ds = A structure containing the parsed dataset.

function ds = parse_dataset(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('file',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addRequired('version',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addOptional('date_format_base','dd/mm/yyyy',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addOptional('date_format_balance','QQ yyyy',@(x)validateattributes(x,{'char'},{'nonempty' 'size' [1 NaN]}));
        ip.addOptional('shares_type','P',@(x)any(validatestring(x,{'P' 'R'})));
        ip.addOptional('crises_type','R',@(x)any(validatestring(x,{'E' 'R'})));
        ip.addOptional('distress_threshold',0.05,@(x)validateattributes(x,{'double'},{'real' 'finite' '>=' 0.05 '<=' 0.75 'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    file = ipr.file;
    version = validate_version(ipr.version);
    date_format_base = validate_date_format(ipr.date_format_base,true);
    date_format_balance = validate_date_format(ipr.date_format_balance,false);
    shares_type = ipr.shares_type;
    crises_type = ipr.crises_type;
    distress_threshold = ipr.distress_threshold;

    file_sheets = validate_xls(ipr.file,'D');

    nargoutchk(1,1);

    ds = parse_dataset_internal(file,file_sheets,version,date_format_base,date_format_balance,shares_type,crises_type,distress_threshold);

end

function ds = parse_dataset_internal(file,file_sheets,version,date_format_base,date_format_balance,shares_type,crises_type,distress_threshold)

    [~,file_name,file_ext] = fileparts(file);
    file_name = [file_name file_ext];

    if (~strcmp(file_sheets{1},'Shares'))
        error(['Error in dataset ''' file_name ''': the first sheet must be the ''Shares'' one.']);
    end

    file_sheets_other = file_sheets(2:end);

    using_prices = strcmp(shares_type,'P');

    if (using_prices)
        tab_shares = parse_table_standard(file,file_name,1,'Shares',date_format_base,[],[],true);
    else
        tab_shares = parse_table_standard(file,file_name,1,'Shares',date_format_base,[],[],false);
    end

    if (any(any(ismissing(tab_shares))))
        error(['Error in dataset ''' file_name ''': the ''Shares'' sheet contains invalid or missing values.']);
    end

    if (width(tab_shares) < 5)
        error(['Error in dataset ''' file_name ''': the ''Shares'' sheet must contain at least the following elements: the observations dates and the time series of the benchmark and at least 3 firms.']);
    end

    n = width(tab_shares) - 2;

    if (using_prices)
        t = height(tab_shares);

        if (t < 253)
            error(['Error in dataset ''' file_name ''': the ''Shares'' sheet must contain at least 253 observations (a full business year plus an additional observation at the beginning of the time series) in order to run consistent calculations.']);
        end

        t = t - 1;
    else
        t = height(tab_shares);

        if (t < 252)
            error(['Error in dataset ''' file_name ''': the ''Shares'' sheet must contain at least 252 observations (a full business year) in order to run consistent calculations.']);
        end
    end

    dates_num = datenum(tab_shares{:,1}); %#ok<DATNM>
    tab_shares.Date = [];

    if (using_prices)
        tab_shares_headers = strtrim(tab_shares.Properties.VariableNames);

        prices = table2array(tab_shares);
        dlp = diff(log(prices));
        dlp(~isfinite(dlp)) = 0;
        prices = prices(:,2:end);

        index_name = tab_shares_headers{1};
        firm_names = tab_shares_headers(2:end);

        index = dlp(:,1);
        returns = dlp(:,2:end);
    else
        tab_shares_headers = strtrim(tab_shares.Properties.VariableNames);

        prices = [];

        index_name = tab_shares_headers{1};
        firm_names = tab_shares_headers(2:end);

        index = tab_shares{:,1};
        returns = tab_shares{:,2:end};
    end

    if (any(cellfun(@isempty,regexpi(firm_names,'^[A-Z][A-Z0-9]{0,11}$'))))
        error(['Error in dataset ''' file_name ''': the ''Shares'' sheet contains invalid firm names (containing invalid characters, not starting with a letter and/or greater than 12 characters).']);
    end

    volumes = [];
    capitalizations = [];

    risk_free_rate = [];
    cds = [];

    assets = [];
    equity = [];
    separate_accounts = [];

    state_variables = [];
    state_variables_names = [];

    groups = 0;
    group_delimiters = [];
    group_names = [];
    group_short_names = [];

    crises = 0;
    crises_dummy = [];
    crisis_dates = [];
    crisis_dummies = [];
    crisis_names = [];

    for tab = {'Volumes' 'Capitalizations' 'CDS' 'Assets' 'Equity' 'Separate Accounts' 'State Variables' 'Groups' 'Crises'}

        tab_index = find(strcmp(file_sheets_other,tab),1);

        if (isempty(tab_index))
            continue;
        end

        tab_index = tab_index + 1;
        tab_name = char(tab);

        switch (tab_name)

            case 'Volumes'

                tab_volumes = parse_table_standard(file,file_name,tab_index,tab_name,date_format_base,dates_num,firm_names,true);
                volumes = table2array(tab_volumes);

            case 'Capitalizations'

                tab_capitalizations = parse_table_standard(file,file_name,tab_index,tab_name,date_format_base,dates_num,firm_names,true);
                capitalizations = table2array(tab_capitalizations);

            case 'CDS'

                tab_cds = parse_table_standard(file,file_name,tab_index,tab_name,date_format_base,dates_num,[{'RF'} firm_names],true);
                tab_arr = table2array(tab_cds);
                risk_free_rate = tab_arr(:,1);
                cds = tab_arr(:,2:end) ./ 10000;

            case 'Assets'

                tab_assets = parse_table_balance(file,file_name,tab_index,tab_name,date_format_balance,dates_num,firm_names,true);
                assets = table2array(tab_assets);

            case 'Equity'

                tab_equity = parse_table_balance(file,file_name,tab_index,tab_name,date_format_balance,dates_num,firm_names,false);
                equity = table2array(tab_equity);

            case 'Separate Accounts'

                tab_separate_accounts = parse_table_balance(file,file_name,tab_index,tab_name,date_format_balance,dates_num,firm_names,true);
                separate_accounts = table2array(tab_separate_accounts);

            case 'State Variables'

                tab_state_variables = parse_table_standard(file,file_name,tab_index,tab_name,date_format_base,dates_num,[],false);
                state_variables = table2array(tab_state_variables);
                state_variables_names = tab_state_variables.Properties.VariableNames;

             case 'Groups'

                [tab_groups,group_counts] = parse_table_groups(file,file_name,tab_index,tab_name,firm_names);
                groups = height(tab_groups);
                group_delimiters = cumsum(group_counts(1:end-1,:));
                group_names = tab_groups{:,1};
                group_short_names = tab_groups{:,2};

             case 'Crises'

                if (strcmp(crises_type,'E'))
                    [tab_crises,crises_dummy] = parse_table_crises_dates(file,file_name,tab_index,tab_name,date_format_base,dates_num);
                    crises = height(tab_crises);
                    crisis_dates = datenum(tab_crises{:,1}); %#ok<DATNM>
                    crisis_names = tab_crises{:,2};
                else
                    [tab_crises,crisis_dummies,crises_dummy] = parse_table_crises_ranges(file,file_name,tab_index,tab_name,date_format_base,dates_num);
                    crises = height(tab_crises);
                    crisis_dates = [datenum(tab_crises{:,2}) datenum(tab_crises{:,3})]; %#ok<DATNM>
                    crisis_names = tab_crises{:,1};
                end

        end
    end

    if (using_prices)
        dates_num = remove_first_observation(dates_num);
        prices = remove_first_observation(prices);
        volumes = remove_first_observation(volumes);
        capitalizations = remove_first_observation(capitalizations);
        risk_free_rate = remove_first_observation(risk_free_rate);
        cds = remove_first_observation(cds);
        assets = remove_first_observation(assets);
        equity = remove_first_observation(equity);
        separate_accounts = remove_first_observation(separate_accounts);
        state_variables = remove_first_observation(state_variables);
        crises_dummy = remove_first_observation(crises_dummy);
        crisis_dummies = remove_first_observation(crisis_dummies);
    end

    [dates_year,~,~,~,~,~] = datevec(dates_num);

    if (~isempty(assets) && ~isempty(equity))
        liabilities = assets - equity;
    else
        liabilities = [];
    end

    if (strcmp(crises_type,'R') && ~isempty(crises_dummy))
        crisis_dates(:,1) = max(crisis_dates(:,1),dates_num(1));
        crisis_dates(:,2) = min(crisis_dates(:,2),dates_num(end));
    end

    if (isempty(prices))
        [defaults,insolvencies] = detect_distress(distress_threshold,returns,volumes,capitalizations,cds,equity);
    else
        [defaults,insolvencies] = detect_distress(distress_threshold,prices,volumes,capitalizations,cds,equity);
    end

    dd_limit = ceil(n * 0.25);
    defaults_cnt = sum(~isnan(defaults));
    insolvencies_cnt = sum(~isnan(insolvencies));

    if (any(defaults < 63))
        error(['Error in dataset ''' file_name ''': it contains firms that defaulted too early and must be removed.']);
    end

    if (defaults_cnt > dd_limit)
        error(['Error in dataset ''' file_name ''': it contains too many defaulted firms (' num2str(defaults_cnt) ' out of a maximum of ' num2str(dd_limit) ').']);
    end

    if (any(insolvencies < 63))
        error(['Error in dataset ''' file_name ''': it contains firms that became insolvent too early and must be removed.']);
    end

    if (insolvencies_cnt > dd_limit)
        error(['Error in dataset ''' file_name ''': it contains too many insolvent firms (' num2str(insolvencies_cnt) ' out of a maximum of ' num2str(dd_limit) ').']);
    end

    includes_volumes = ismember({'Volumes'},file_sheets_other);
    includes_capitalizations = ismember({'Capitalizations'},file_sheets_other);
    includes_cds = ismember({'CDS'},file_sheets_other);
    includes_balance_sheet = sum(ismember({'Assets' 'Equity'},file_sheets_other)) == 2;
    includes_crises = ~isempty(crises_dummy);

    if (using_prices && includes_capitalizations)
        supports_bubbles_detection = true;
    else
        supports_bubbles_detection = false;
        warning('MATLAB:SystemicRisk',['The dataset ''' file_name ''' does not meet all the requirements for the computation of bubbles detection measures (''Shares'' must be expressed as prices, the ''Capitalizations'' sheet must be included).']);
    end

    if (includes_cds)
        supports_cross_entropy = true;
    else
        supports_cross_entropy = false;
        warning('MATLAB:SystemicRisk',['The dataset ''' file_name ''' does not meet all the requirements for the computation of cross-entropy measures (the ''CDS'' sheet must be included).']);
    end

    if (includes_capitalizations && includes_balance_sheet)
        supports_cross_sectional = true;
    else
        supports_cross_sectional = false;
        warning('MATLAB:SystemicRisk',['The dataset ''' file_name ''' does not meet all the requirements for the computation of cross-sectional measures (''Assets'', ''Capitalizations'' and ''Equity'' sheets must be included).']);
    end

    if (includes_capitalizations && includes_cds && includes_balance_sheet)
        supports_default = true;
    else
        supports_default = false;
        warning('MATLAB:SystemicRisk',['The dataset ''' file_name ''' does not meet all the requirements for the computation of default measures (''Assets'', ''Capitalizations'', ''CDS'' and ''Equity'' sheets must be included).']);
    end

    if (using_prices && includes_volumes && includes_capitalizations)
        supports_liquidity = true;
    else
        supports_liquidity = false;
        warning('MATLAB:SystemicRisk',['The dataset ''' file_name ''' does not meet all the requirements for the computation of liquidity measures (''Shares'' must be expressed as prices, ''Capitalizations'' and ''Volumes'' sheets must be included).']);
    end

    if (includes_crises)
        supports_comparison = true;
    else
        supports_comparison = false;
        warning('MATLAB:SystemicRisk',['The dataset ''' file_name ''' does not meet all the requirements for the comparison of systemic risk measures (the ''Crises'' sheet must be included).']);
    end

    returns = distress_data(returns,defaults);
    prices = distress_data(prices,defaults);
    volumes = distress_data(volumes,defaults);
    capitalizations = distress_data(capitalizations,defaults);
    cds = distress_data(cds,defaults);
    assets = distress_data(assets,defaults);
    equity = distress_data(equity,defaults);
    liabilities = distress_data(liabilities,defaults);
    separate_accounts = distress_data(separate_accounts,defaults);

    ds = struct();

    ds.File = file;
    ds.Version = version;
    ds.CreationDate = now(); %#ok<TNOW1>

    ds.Result = [];
    ds.ResultDate = [];
    ds.ResultAnalysis = [];
    ds.ResultSerial = [];

    ds.TimeSeries = {'Assets' 'Capitalizations' 'CDS' 'Equity' 'Liabilities' 'Returns' 'SeparateAccounts' 'Volumes'};

    ds.N = n;
    ds.T = t;

    ds.DatesNum = dates_num;
    ds.DatesStr = cellstr(datestr(datetime(dates_num,'ConvertFrom','datenum'),date_format_base)); %#ok<DATST>
    ds.MonthlyTicks = numel(unique(dates_year)) <= 3;

    ds.IndexName = index_name;
    ds.FirmNames = firm_names;

    ds.Index = index;
    ds.Returns = returns;

    ds.Prices = prices;
    ds.Volumes = volumes;
    ds.Capitalizations = capitalizations;

    ds.RiskFreeRate = risk_free_rate;
    ds.CDS = cds;

    ds.Assets = assets;
    ds.Equity = equity;
    ds.Liabilities = liabilities;
    ds.SeparateAccounts = separate_accounts;

    ds.StateVariables = state_variables;
    ds.StateVariablesNames = state_variables_names;

    ds.Groups = groups;
    ds.GroupDelimiters = group_delimiters;
    ds.GroupNames = group_names;
    ds.GroupShortNames = group_short_names;

    ds.Crises = crises;
    ds.CrisesType = crises_type;
    ds.CrisesDummy = crises_dummy;
    ds.CrisisDates = crisis_dates;
    ds.CrisisDummies = crisis_dummies;
    ds.CrisisNames = crisis_names;

    ds.Defaults = defaults;
    ds.Insolvencies = insolvencies;

    ds.SupportsBubblesDetection = supports_bubbles_detection;
    ds.SupportsComponent = true;
    ds.SupportsConnectedness = true;
    ds.SupportsCrossEntropy = supports_cross_entropy;
    ds.SupportsCrossQuantilogram = true;
    ds.SupportsCrossSectional = supports_cross_sectional;
    ds.SupportsDefault = supports_default;
    ds.SupportsLiquidity = supports_liquidity;
    ds.SupportsRegimeSwitching = true;
    ds.SupportsSpillover = true;
    ds.SupportsTailDependence = true;
    ds.SupportsComparison = supports_comparison;

end

%% VALIDATION

function date_format = validate_date_format(date_format,base)

    try
        datestr(now(),date_format); %#ok<DATST,TNOW1>
    catch e
        error(['The date format ''' date_format ''' is invalid.' new_line() strtrim(regexprep(e.message,'Format.+$',''))]);
    end

    if (~any(regexp(date_format,'(?<!y)(y{2}|y{4})(?!y)')))
        error('The date format must always include year information.');
    end

    if (any(regexp(date_format,'(?:HH|MM|SS|FFF|AM|PM)')))
        error('The date format must not include time information.');
    end

    if (base && ~any(regexp(date_format,'(?<!d)(d{2}|d{4})(?!d)')))
        error('The base date format must be tied to a daily frequency.');
    end

end

function version = validate_version(version)

    if (~any(regexp(version,'v\d+\.\d+')))
        error('The version format is invalid.');
    end

end

%% PARSING

function tab = ensure_field_consistency(file_name,name,tab,index,output_type,target_type,date_format_dt)

    switch (target_type)

        case 'char'

            if (strcmp(output_type,'cell'))
                field = vertcat(tab{:,index});

                if (~ischar(field{1}))
                    error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains invalid column types.']);
                end

                tab.(tab.Properties.VariableNames{index}) = field;
            elseif (strcmp(output_type,'char'))
                tab.(tab.Properties.VariableNames{index}) = cellstr(tab{:,index});
            else
                error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains invalid column types.']);
            end

        case 'datetime'

            if (strcmp(output_type,'cell'))
                if (iscell(tab{1,index}))
                    field = cellfun(@(x)x,tab{:,index},'UniformOutput',false);
                end

                if (ischar(field{1}))
                    field = cellfun(@(x)datetime(x,'InputFormat',date_format_dt),field,'UniformOutput',false);
                elseif (~isdatetime(field{1}))
                    error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains invalid column types.']);
                end

                tab.(tab.Properties.VariableNames{index}) = [field{:}].';
            elseif (strcmp(output_type,'char'))
                field = cellfun(@(x)datetime(x,'InputFormat',date_format_dt),cellstr(tab{:,index}));
                tab.(tab.Properties.VariableNames{1}) = field;
            elseif (~strcmp(output_type,'datetime'))
                error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains invalid column types.']);
            end

        case 'double'

            if (strcmp(output_type,'cell'))
                if (iscell(tab{1,index}))
                    field = cellfun(@(x)x{:},tab{:,index},'UniformOutput',false);
                end

                field_first = field{1};

                if (~isa(field_first,'double') || ~isscalar(field_first))
                    error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains invalid column types.']);
                end

                tab.(tab.Properties.VariableNames{index}) = [field{:}].';
            elseif (~strcmp(output_type,'double'))
                error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains invalid column types.']);
            end

        otherwise

            error(['Invalid target type specified: ' target_type '.']);

    end

end

function tab = parse_table_balance(file,file_name,index,name,date_format,dates_num,firm_names,check_negatives)

    data_types = [{'datetime'} repmat({'double'},1,numel(firm_names))];
    tab_partial = read_table(file,file_name,index,name,date_format,data_types);

    if (~strcmp(tab_partial.Properties.VariableNames(1),'Date'))
        error(['Error in dataset ''' file_name ''': the first column of the ''' name ''' sheet must be called ''Date'' and must contain the observation dates.']);
    end

    if (any(any(~isfinite(tab_partial{:,2:end}))))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains invalid values.']);
    end

    if (check_negatives && any(any(tab_partial{:,2:end} < 0)))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains negative values.']);
    end

    t_current = height(tab_partial);
    dates_num_current = datenum(tab_partial.Date); %#ok<DATNM>
    tab_partial.Date = [];

    if (t_current ~= numel(unique(dates_num_current)))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains duplicate observation dates.']);
    end

    if (any(dates_num_current ~= sort(dates_num_current)))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains unsorted observation dates.']);
    end

    if (t_current ~= numel(unique(cellstr(datestr(dates_num,date_format))))) %#ok<DATST>
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains an invalid number of observation dates.']);
    end

    if (~isequal(tab_partial.Properties.VariableNames,firm_names))
        error(['Error in dataset ''' file_name ''': the firm names between the ''Shares'' sheet and the ''' name ''' sheet are mismatching.']);
    end

    dates_from = cellstr(datestr(dates_num_current,date_format)).'; %#ok<DATST>
    dates_to = cellstr(datestr(dates_num,date_format)).'; %#ok<DATST>

    if (any(~ismember(dates_to,dates_from)))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet observation dates do not cover all the ''Shares'' sheet observation dates.']);
    end

    t = numel(dates_num);
    n = numel(firm_names);
    tab = array2table(NaN(t,n),'VariableNames',tab_partial.Properties.VariableNames);

    for date = dates_from
        members_from = ismember(dates_from,date);
        members_to = ismember(dates_to,date);
        tab{members_to,:} = repmat(tab_partial{members_from,:},sum(members_to),1);
    end

end

function [tab,crises_dummy] = parse_table_crises_dates(file,file_name,index,name,date_format,dates_num)

    data_types = {'datetime' 'char'};
    tab = read_table(file,file_name,index,name,date_format,data_types);

    if (~isequal(tab.Properties.VariableNames,{'Date' 'Name'}))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains invalid (wrong name) or misplaced (wrong order) columns.']);
    end

    dates_num_min = min(dates_num);
    dates_num_max = max(dates_num);

    dates_num_current = datenum(tab.Date); %#ok<DATNM>

    if (any(dates_num_current < dates_num_min) || any(dates_num_current > dates_num_max))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains crisis events outside the scope of the dataset (' datestr(dates_num_min,date_format) ' - ' datestr(dates_num_max,date_format) ').']); %#ok<DATST>
    end

    n = height(tab);
    crises_dummy = zeros(numel(dates_num),1);

    for i = 1:n
        [check,index] = ismember(datenum(tab.Date(i)),dates_num); %#ok<DATNM>

        if (~check)
            error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains crisis events outside the scope of the dataset.']);
        end

        crises_dummy(index) = 1;
    end

    if (all(crises_dummy == 1))
        error(['Error in dataset ''' file_name ''': all the observation are considered to be part of a distress period.']);
    end

end

function [tab,crisis_dummies,crises_dummy] = parse_table_crises_ranges(file,file_name,index,name,date_format,dates_num)

    data_types = {'char' 'datetime' 'datetime'};
    tab = read_table(file,file_name,index,name,date_format,data_types);

    if (~isequal(tab.Properties.VariableNames,{'Name' 'StartDate' 'EndDate'}))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains invalid (wrong name) or misplaced (wrong order) columns.']);
    end

    if (any(cellfun(@length,tab.Name) > 30))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains one or more crises with a name longer than 25 characters.']);
    end

    dates_num_min = min(dates_num);
    dates_num_max = max(dates_num);

    dates_num_start = datenum(tab.StartDate); %#ok<DATNM>
    dates_num_end = datenum(tab.EndDate); %#ok<DATNM>

    if (any(dates_num_start > dates_num_max) || any(dates_num_end < dates_num_min))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains crises outside the scope of the dataset (' datestr(dates_num_min,date_format) ' - ' datestr(dates_num_max,date_format) ').']); %#ok<DATST>
    end

    if (any(dates_num_start > dates_num_end))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains one or more crises with start date greater than end date.']);
    end

    n = height(tab);
    crisis_dummies = zeros(numel(dates_num),n);

    for i = 1:n
        seq = datenum(tab.StartDate(i):tab.EndDate(i)).'; %#ok<DATNM>
        crisis_dummies(ismember(dates_num,seq),i) = 1;
    end

    crises_dummy = double(sum(crisis_dummies,2) > 0);

    if (all(crises_dummy == 1))
        error(['Error in dataset ''' file_name ''': all the observation are considered to be part of a distress period.']);
    end

end

function [tab,group_counts] = parse_table_groups(file,file_name,index,name,firm_names)

    data_types = {'char' 'char' 'double'};
    tab = read_table(file,file_name,index,name,[],data_types);

    if (~isequal(tab.Properties.VariableNames,{'Name' 'ShortName' 'Count'}))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains invalid (wrong name) or misplaced (wrong order) columns.']);
    end

    if (any(~isfinite(tab{:,end})))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains invalid values.']);
    end

    group_counts = tab.Count;

    if (any(group_counts <= 0) || any(round(group_counts) ~= group_counts))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains one or more groups with an invalid number of firms.']);
    end

    if (sum(group_counts) ~= numel(firm_names))
        error(['Error in dataset ''' file_name ''': in the ''' name ''' sheet, the number of firms must be equal to the one defined in the ''Shares'' sheet.']);
    end

    tab.Name = strtrim(tab.Name);
    tab.ShortName = strtrim(tab.ShortName);

    if (any(cellfun(@isempty,regexpi(tab.ShortName,'^[A-Z][A-Z0-9]{0,4}$'))))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains groups with an invalid short name (containing invalid characters, not starting with a letter and/or greater than 5 characters).']);
    end

end

function tab = parse_table_standard(file,file_name,index,name,date_format,dates_num,firm_names,check_negatives)

    if (isempty(firm_names))
        data_types = [];
    else
        data_types = [{'datetime'} repmat({'double'},1,numel(firm_names))];
    end

    tab = read_table(file,file_name,index,name,date_format,data_types);

    if (~strcmp(tab.Properties.VariableNames(1),'Date'))
        error(['Error in dataset ''' file_name ''': the first column of the ''' name ''' sheet must be called ''Date'' and must contain the observation dates.']);
    end

    if (any(any(~isfinite(tab{:,2:end}))))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains invalid values.']);
    end

    if (check_negatives)
        if (~isempty(firm_names) && strcmp(firm_names{1},'RF'))
            if (any(any(tab{:,3:end} < 0)))
                error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains negative values.']);
            end
        else
            if (any(any(tab{:,2:end} < 0)))
                error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains negative values.']);
            end
        end
    end

    t_current = height(tab);
    dates_num_current = datenum(tab.Date); %#ok<DATNM>

    if (t_current ~= numel(unique(dates_num_current)))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains duplicate observation dates.']);
    end

    if (any(dates_num_current ~= sort(dates_num_current)))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains unsorted observation dates.']);
    end

    if (~isempty(dates_num))
        if ((t_current ~= numel(dates_num)) || any(dates_num_current ~= dates_num))
            error(['Error in dataset ''' file_name ''': the observation dates between the ''Shares'' sheet and the ''' name ''' sheet are mismatching.']);
        end

        tab.Date = [];
    end

    if (~isempty(firm_names) && ~isequal(tab.Properties.VariableNames,firm_names))
        error(['Error in dataset ''' file_name ''': the firm names between the ''Shares'' sheet and the ''' name ''' sheet are mismatching.']);
    end

end

function tab = read_table(file,file_name,index,name,date_format,data_types)

    date_format_dt = date_format;
    date_format_dt = strrep(date_format_dt,'m','M');
    date_format_dt = strrep(date_format_dt,'QQ','QQQ');

    data_types_count = numel(data_types);

    if (verLessThan('MATLAB','9.1'))
        try
            if (ispc())
                try
                    tab = readtable(file,'Sheet',index,'Basic',true);
                catch
                    tab = readtable(file,'Sheet',index);
                end
            else
                tab = readtable(file,'Sheet',index);
            end
        catch e
            error(['Error in dataset ''' file_name ''': the ''' name ''' sheet could not be read.' new_line() e.message]);
        end

        output_vars = varfun(@class,tab,'OutputFormat','cell');

        if (data_types_count == 0)
            data_types_count = width(tab);

            if (data_types_count > 0)
                data_types = [{'datetime'} repmat({'double'},1,data_types_count - 1)];
            end
        end

        for i = 1:data_types_count
            tab = ensure_field_consistency(file_name,name,tab,i,output_vars{i},data_types{i},date_format_dt);
        end
    else
        if (ispc())
            options = detectImportOptions(file,'Sheet',index,'NumHeaderLines',0);
        else
            options = detectImportOptions(file,'Sheet',name,'NumHeaderLines',0);
        end

        if (data_types_count == 0)
            data_types_count = numel(options.VariableNames);

            if (data_types_count > 0)
                data_types = [{'datetime'} repmat({'double'},1,data_types_count - 1)];
            end
        end

        options = setvartype(options,options.VariableNames(1:data_types_count),data_types);

        dt_indices = strcmp(data_types,'datetime');

        if (any(dt_indices))
            dt_variables = options.VariableNames(dt_indices);
            options = setvaropts(options,dt_variables,'InputFormat',date_format_dt);
        end

        try
            if (ispc())
                try
                    tab = readtable(file,options,'Basic',true);
                catch
                    tab = readtable(file,options);
                end
            else
                tab = readtable(file,options);
            end
        catch e
            error(['Error in dataset ''' file_name ''': the ''' name ''' sheet could not be read.' new_line() e.message]);
        end
    end

    h = height(tab);
    w = width(tab);

    if (h == 0)
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet must contain at least 1 row.']);
    end

    imr = all(ismissing(tab),2);
    f = find(diff([false; imr; false] == 0));

    if (~isempty(f))
        indices = f(1:2:end-1);
        counts = f(2:2:end) - indices;

        index_last = indices(end);
        count_last = counts(end);

        if ((index_last + count_last - 1) == h)
            tab(index_last:end,:) = [];
        end
    end

    isc = all(ismissing(tab),1).';
    f = find(diff([false; isc; false] == 0));

    if (~isempty(f))
        indices = f(1:2:end-1);
        counts = f(2:2:end) - indices;

        index_last = indices(end);
        count_last = counts(end);

        if ((index_last + count_last - 1) == w)
            tab(:,index_last:end) = [];
        end
    end

    if ((data_types_count > 1) && (width(tab) ~= data_types_count))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet must contain exactly ' num2str(data_types_count) ' columns.']);
    end

    if (~all(cellfun(@isempty,regexp(tab.Properties.VariableNames,'^Var\d+$','once'))))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains unnamed columns.']);
    end

    if (any(any(ismissing(tab))))
        error(['Error in dataset ''' file_name ''': the ''' name ''' sheet contains missing values.']);
    end

end

%% COMPUTATIONS

function [defaults,insolvencies] = detect_distress(distress_threshold,shares,volumes,capitalizations,cds,equity)

    [t,n] = size(shares);
    limit = round(t * distress_threshold,0);

    data = {shares volumes capitalizations cds};
    k = sum(~cellfun(@isempty,data,'UniformOutput',true));

    data = [shares volumes capitalizations cds];
    defaults = NaN(1,n);

    for i = 1:n
        data_i = data(:,i:n:(n * k));
        defaults_i = NaN(1,k);

        for j = 1:k
            data_ij = data_i(:,j);

            f = find(diff([1; data_ij; 1] == 0));

            if (~isempty(f))
                indices = f(1:2:end-1);
                counts = f(2:2:end) - indices;

                index_last = indices(end);
                count_last = counts(end);

                if (((index_last + count_last - 1) == t) && (count_last >= limit))
                    defaults_i(j) = index_last;
                end
            end
        end

        defaults(i) = min(defaults_i,[],'omitnan');
    end

    insolvencies = NaN(1,n);

    if (~isempty(equity))
        for i = 1:n
            eq = equity(:,i);

            f = find(diff([false; eq < 0; false] ~= 0));

            if (~isempty(f))
                indices = f(1:2:end-1);
                counts = f(2:2:end) - indices;

                index_last = indices(end);
                count_last = counts(end);

                if (((index_last + count_last - 1) == numel(eq)) && (count_last >= limit))
                    insolvencies(i) = index_last;
                end
            end
        end
    end

end

function data = remove_first_observation(data)

    if (~isempty(data))
        data = data(2:end,:);
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
classdef NonLinear_Solver < handle

    properties (Access = public)
    end

    properties (Access = private)
        tol
        state
        solver
        element
        free_dof
    end

    methods (Access = public)

        function obj = NonLinear_Solver(cParams)
            obj.init(cParams);
            obj.createSolver();
        end

        function x = solve(obj, params)
            switch params.state
                case 'Steady'
                    x = obj.solveSteady();
                case 'Transient'
                    x = obj.solveTransient(params);
            end
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.tol      = cParams.tol;
%             obj.state    = cParams.state;
            obj.element  = cParams.element;
            obj.free_dof = cParams.free_dof;
        end

        function createSolver(obj)
            s.type = 'DIRECT';
            obj.solver = Solver.create(s);
        end

        function sol = solveSteady(obj)
            total_free_dof = sum(obj.free_dof);
            dt = Inf;
            dr = obj.element.computedr(dt);
            x0 = zeros(total_free_dof,1);

            r = obj.element.computeResidual(x0,dr);
            x = obj.convergeSolution(dr, r, x0);
            sol = x;
        end

        function sol = solveTransient(obj,params)
            dt = params.dt;
            final_time = params.final_time;
            total_free_dof = sum(obj.free_dof);
            x_n(:,1) = zeros(total_free_dof,1);
            x0 = zeros(total_free_dof,1);

            dr = obj.element.computedr(dt);

            for istep = 2: final_time/dt
                u_previous_step = x_n(1:obj.free_dof(1),istep-1);

                r = obj.element.computeResidual(x0,dr,u_previous_step);
                while dot(r,r) > obj.tol
                    inc_x = obj.solver.solve(dr,-r);
                    x = x0 + inc_x;
                    % Compute r
                    r = obj.element.computeResidual(x,dr,u_previous_step);
                    x0 = x;
                end
                x_n(:,istep) = x;
            end
            sol = x_n;
        end

        function sol = convergeSolution(obj, dr, r, x0, u_previous_step)
            while dot(r,r) > obj.tol
                inc_x = obj.solver.solve(dr,-r);
                x = x0 + inc_x;
                % Compute r
                r = obj.element.computeResidual(x,dr);
                x0 = x;
            end
            sol = x0;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function han = m_scatter(lon, lat, varargin)
% function m_scatter(lon, lat, S, C)

m_proj('miller','lat',82);
%m_grid('linestyle','none','box','on','tickdir','out');
[x,y] = m_ll2xy(lon, lat);
%drawnow
h = scatter(x, y, varargin{:});
line(x(1:10), y(1:10))
%plot(x,y,'.')
%get(gca,
z = findobj(gca, 'type', 'scattergroup')
m_grid('linestyle','none','box','on','tickdir','out')
%h = scatter(x, y, varargin{:});

if nargout == 1
  han = h;
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
classdef testMaterialInsertingLpBall < testShowingError

    properties (Access = protected)
        tol = 5e-2;
        testName = 'testMaterialInsertingLpBall';
    end

    properties (Access = private)
        vadSAcomparator
        plotter
    end

    methods (Access = public)

        function obj = testMaterialInsertingLpBall()
            obj.computeVademecumAndSimpAllValues();
            obj.plotVademecumAndSimpAllValues();
        end

    end

    methods (Access = protected)

        function computeError(obj)
            obj.error = 0;
        end

    end

    methods (Access = private)

        function computeVademecumAndSimpAllValues(obj)
            s.fileName =  'VademecumSmoothCorner';
            v = VademecumSimpAllComparator(s);
            v.calculate();
            obj.vadSAcomparator = v;
        end

        function plotVademecumAndSimpAllValues(obj)
            s.vadSAcomparator = obj.vadSAcomparator;
            obj.plotter = VademecumComparatorPlotter(s);
            obj.plotter.plot();
        end

    end

end


text:
%DTRANSFORM Distance transform
%
% DT = DTRANSFORM(IM, OPTIONS) is the distance transform of the
% binary image IM. The value of each output pixel is the distance (pixels)
% to the closest set pixel.
%
% Options::
% 'Euclidean'   use Euclidean distance (default)
% 'cityblock'   use cityblock (Manhattan) distance
% 'show',T      display the evolving distance transform, with a delay of T
%               seconds between frames
%
% See also IMORPH, DISTANCEXFORM, DXform.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function d = dtransform(world, varargin)

    opt.metric = {'Euclidean', 'cityblock'};
    opt.show = [];

    [opt,args] = tb_optparse(opt, varargin);
    if length(args) > 0 && isnumeric(args{1})
        opt.show = args{1};
    end

    if strcmpi(opt.metric, 'cityblock')
        m = ones(3,3);
        m(2,2) = 0;
    elseif strcmpi(opt.metric, 'Euclidean')
        r2 = sqrt(2);
        m = [r2 1 r2; 1 0 1; r2 1 r2];
    end

    world(world==0) = Inf;
    world(world==1) = 0;

    count = 0;
    while 1
        world = imorph(world, m, 'plusmin');
        count = count+1;
        if opt.show
            cmap = gray(256);
            cmap = [1 0 0; cmap];
            colormap(cmap)
            image(world+1, 'CDataMapping', 'direct');
            set(gca, 'Ydir', 'normal');
            xlabel('x');
            ylabel('y');
            pause(opt.show);
        end

        if length(find(world(:)==Inf)) == 0
            break;
        end
    end

    if opt.show
        fprintf('%d iterations\n', count);
    end

    d = world;

author: petercorke

repo: machinevision-toolbox-matlab


text:
function []=ps_sb_merge()
%PS_SB_MERGE merge PS and SB pixels
%
%   Andy Hooper May 2007
%
logit;

ps_wd=pwd;
i=strfind(ps_wd,'INSAR_');
if isempty(i)
    error('INSAR_* does not appear in current directory path')
end

sb_wd1=[ps_wd(1:i+13),'/SMALL_BASELINES'];
sb_wd=[sb_wd1,ps_wd(i+14:end)];
merged_wd1=[ps_wd(1:i+13),'/MERGED'];
merged_wd=[merged_wd1,ps_wd(i+14:end)];

if ~exist(merged_wd)
    mkdir(merged_wd)
end
if ~exist([merged_wd,'/patch_noover.in']) & exist([sb_wd,'/patch_noover.in'])
    aa=['!cp ',sb_wd,'/patch_noover.in ',merged_wd];
    eval(aa)
end
if ~exist([merged_wd1,'/parms.mat'])
    aa=['!cp ',sb_wd1,'/parms.mat ',merged_wd1,'/parms.mat'];
    eval(aa)
end

psver=2;
save([merged_wd,'/psver'],'psver')

psps=load([ps_wd,'/ps2.mat']);
pssb=load([sb_wd,'/ps2.mat']);

pssb_ix=zeros(psps.n_ps,1); % points from ps to sb

imax=max([psps.ij(:,2);pssb.ij(:,2)])+1;
jmax=max([psps.ij(:,3);pssb.ij(:,3)])+1;
grid_ix=zeros(imax,jmax);
grid_ix(pssb.ij(:,2)+1+pssb.ij(:,3)*imax)=[1:pssb.n_ps];
pssb_ix=grid_ix(psps.ij(:,2)+1+psps.ij(:,3)*imax);

psu_ix=find(pssb_ix==0); % ps only index
psnu_ix=find(pssb_ix);  % ps non-unique index
sbnu_ix=pssb_ix(psnu_ix); % sb non-unique index

%%%%%%%%% Some non-adjacent pixels are allocated the same lon/lat by DORIS. %%%%%%%%%
[dups,I]=setdiff(psps.lonlat(psu_ix,:),pssb.lonlat,'rows');
fprintf('%d pixels with duplicate lon/lat dropped\n',length(psu_ix)-length(dups))
psu_ix=psu_ix(I);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

G=zeros(pssb.n_ifg,pssb.n_image);
for i=1:pssb.n_ifg
     G(i,pssb.ifgday_ix(i,1))=-1;
     G(i,pssb.ifgday_ix(i,2))=1;
end

pmps=load([ps_wd,'/pm2.mat']);
pmsb=load([sb_wd,'/pm2.mat']);

pmps.ph_res_all=[pmps.ph_res(:,1:psps.master_ix-1),pmps.C_ps,pmps.ph_res(:,psps.master_ix:end)]; % add master residual
pmps.ph_res3=exp(j*(G*pmps.ph_res_all')'); % would-be residuals of sb ifgs
pmps.coh_ps2=abs(sum(pmps.ph_res3,2))/pssb.n_ifg; % would-be coherence of sb ifgs

psnu_coh=pmps.coh_ps2(psnu_ix);
psnu_snr=1./(1./psnu_coh.^2-1);
sbnu_coh=pmsb.coh_ps(sbnu_ix);
sbnu_snr=1./(1./sbnu_coh.^2-1);
ps_high_coh_ix=pmps.coh_ps2(psu_ix)>min(pmsb.coh_ps);
psu_ix=psu_ix(ps_high_coh_ix);
psu_coh=pmps.coh_ps2(psu_ix);

save([merged_wd,'/merge2'],'psu_ix','psu_coh','psnu_ix','sbnu_ix','psnu_coh','sbnu_coh','G')

ps=pssb;
ps.ij=[psps.ij(psu_ix,:);pssb.ij];
ps.lonlat=[psps.lonlat(psu_ix,:);pssb.lonlat];
ll0=(max(ps.lonlat)+min(ps.lonlat))/2;
xy=llh2local(ps.lonlat',ll0)*1000;
xy=xy';
sort_x=sortrows(xy,1);
sort_y=sortrows(xy,2);
n_pc=round(size(xy,1)*0.001);
bl=mean(sort_x(1:n_pc,:)); % bottom left corner
tr=mean(sort_x(end-n_pc:end,:)); % top right corner
br=mean(sort_y(1:n_pc,:)); % bottom right  corner
tl=mean(sort_y(end-n_pc:end,:)); % top left corner

heading=getparm('heading');
theta=(180-heading)*pi/180;
if theta>pi
    theta=theta-2*pi;
end

rotm=[cos(theta),sin(theta); -sin(theta),cos(theta)];
xy=xy';
xynew=rotm*xy; % rotate so that scene axes approx align with x=0 and y=0
if max(xynew(1,:))-min(xynew(1,:))<max(xy(1,:))-min(xy(1,:)) &...
   max(xynew(2,:))-min(xynew(2,:))<max(xy(2,:))-min(xy(2,:))
    xy=xynew; % check that rotation is an improvement
    disp(['Rotating by ',num2str(theta*180/pi),' degrees']);
end

xy=single(xy');
[dummy,sort_ix]=sortrows(xy,[2,1]); % sort in ascending y order
xy=xy(sort_ix,:);
xy=[[1:size(xy,1)]',xy];
xy(:,2:3)=round(xy(:,2:3)*1000)/1000; % round to mm
ps.xy=xy;
ps.lonlat=ps.lonlat(sort_ix,:);
ps.ij=ps.ij(sort_ix,:);
ps.n_ps=size(ps.xy,1);
save([merged_wd,'/ps2'],'-struct','ps')

pmps.ph_patch2=exp(j*G*angle([pmps.ph_patch(:,1:psps.master_ix-1),ones(psps.n_ps,1),pmps.ph_patch(:,psps.master_ix:end)])').';
pm.ph_patch=pmsb.ph_patch;
pm.ph_patch(sbnu_ix,:)=pmsb.ph_patch(sbnu_ix,:).*repmat(sbnu_snr,1,pssb.n_ifg)+pmps.ph_patch2(psnu_ix,:).*repmat(psnu_snr,1,pssb.n_ifg);
pm.ph_patch(sbnu_ix,:)=pm.ph_patch(sbnu_ix,:)./abs(pm.ph_patch(sbnu_ix,:));
pm.ph_patch=[pmps.ph_patch2(psu_ix,:);pm.ph_patch];
pm.ph_patch=pm.ph_patch(sort_ix,:);

save([merged_wd,'/pm2'],'-struct','pm')
clear pm pmps pmsb

rcps=load([ps_wd,'/rc2.mat']);
rcsb=load([sb_wd,'/rc2.mat']);
rcps.ph_rc2=exp(j*G*angle(rcps.ph_rc)').'; % convert phase to sb ifgs
rc=rcsb;
rc.ph_rc(sbnu_ix,:)=rcsb.ph_rc(sbnu_ix,:).*repmat(sbnu_snr,1,pssb.n_ifg)+rcps.ph_rc2(psnu_ix,:).*repmat(psnu_snr,1,pssb.n_ifg);
rc.ph_rc(sbnu_ix,:)=rc.ph_rc(sbnu_ix,:)./abs(rc.ph_rc(sbnu_ix,:));
%rc.ph_rc(sbnu_ix,:)=rc.ph_rc(sbnu_ix,:)+rcps.ph_rc2(psnu_ix,:);
rc.ph_rc=[rcps.ph_rc2(psu_ix,:);rc.ph_rc];
rc.ph_rc=rc.ph_rc(sort_ix,:);

save([merged_wd,'/rc2'],'-struct','rc')
clear rc rcps rcsb

bpps=load([ps_wd,'/bp2.mat']);
bpsb=load([sb_wd,'/bp2.mat']);
bpps.bperp_mat2=(G*[bpps.bperp_mat(:,1:psps.master_ix-1),zeros(psps.n_ps,1),bpps.bperp_mat(:,psps.master_ix:end)]')';
bp=bpsb;
bp.bperp_mat=[bpps.bperp_mat2(psu_ix,:);bp.bperp_mat];
bp.bperp_mat=bp.bperp_mat(sort_ix,:);
save([merged_wd,'/bp2'],'-struct','bp')
clear bp bpps bpsb

hgtps=load([ps_wd,'/hgt2.mat']);
hgtsb=load([sb_wd,'/hgt2.mat']);
hgt=hgtsb;
hgt.hgt=[hgtps.hgt(psu_ix,:);hgt.hgt];
hgt.hgt=hgt.hgt(sort_ix,:);
save([merged_wd,'/hgt2'],'-struct','hgt')
clear hgt hgtps hgtsb

laps=load([ps_wd,'/la2.mat']);
lasb=load([sb_wd,'/la2.mat']);
la=lasb;
la.la=[laps.la(psu_ix,:);la.la];
la.la=la.la(sort_ix,:);
save([merged_wd,'/la2'],'-struct','la')
clear la laps lasb
logit(1);

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function [imo,rois] = fast_rcnn_eval_get_batch(images, imdb, batch, opts)
% FAST_RCNN_GET_BATCH_EVAL Load, preprocess, and pack images for CNN
% evaluation

% opts.numFgRoisPerImg = 128;
% opts.numRoisPerImg = 64;
% opts.maxScale = 1000;
% opts.bgLabel = 21;
% opts.visualize = 1;
% opts.scale = 600;
% opts.interpolation = 'bicubic';
% opts.averageImage = [];
% opts.numThreads = 2;
%
% Copyright (C) 2016 Hakan Bilen.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if isempty(images)
  imo = [] ;
  rois = [] ;
  return ;
end

% fetch is true if images is a list of filenames (instead of
% a cell array of images)
fetch = ischar(images{1}) ;

% prefetch is used to load images in a separate thread
prefetch = fetch & opts.prefetch ;

if prefetch
  vl_imreadjpeg(images, 'numThreads',opts.numThreads,'prefetch') ;
  imo = [] ;
  rois = [] ;
  return ;
end

if fetch
  ims = vl_imreadjpeg(images,'numThreads',opts.numThreads) ;
else
  ims = images ;
end

imre = cell(1,numel(batch));
maxW = 0;
maxH = 0;

pboxes  = cell(1,numel(batch));

% get fg and bg rois
for b=1:numel(batch)
  pbox   = imdb.boxes.pbox{batch(b)};

  if size(pbox,2)~=4
    error('wrong box size');
  end

  pboxes{b} = pbox;
end

% rescale images and rois
rois = [];
for b=1:numel(batch)
  imSize = size(ims{b});

  h = imSize(1);
  w = imSize(2);

  factor = max(opts.scale(1)/h,opts.scale(1)/w);

  if any([h*factor,w*factor]>opts.maxScale)
    factor = min(opts.maxScale/h,opts.maxScale/w);
  end

  if abs(factor-1)>1e-3
    imre{b} = imresize(ims{b},factor,'Method',opts.interpolation,...
        'antialiasing', false);
  else
    imre{b} = ims{b};
  end

  if imdb.boxes.flip(batch(b))
    im = imre{b};
    imre{b} = im(:,end:-1:1,:);
  end

  imreSize = size(imre{b});

  maxH = max(imreSize(1),maxH);
  maxW = max(imreSize(2),maxW);

  % adapt bounding boxes into new coord
  bbox = pboxes{b};
  if any(bbox(:)<=0)
    error('bbox error');
  end

  nB = size(bbox,1);
  tbbox = bbox_scale(bbox,factor,[imreSize(2) imreSize(1)]);
  if any(tbbox(:)<=0)
    error('tbbox error');
  end

  rois = [rois [b*ones(1,nB); tbbox' ] ];
end

% rois = single(rois);
imo = zeros(maxH,maxW,size(imre{1},3),numel(batch),'single');
for b=1:numel(batch)
  if ~isempty(opts.averageImage)
    imre{b} = single(bsxfun(@minus,imre{b},opts.averageImage));
  end
  sz = size(imre{b});
  imo(1:sz(1),1:sz(2),:,b) = single(imre{b});
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function semantic_segmentation_rgb(fileName, sensorType, varName)

    %% Check if input is valid
    if ~isstring(fileName) && ~ischar(fileName)
        error("File name must a string or char array\n")
    end

    if ~strcmp(sensorType, "semantic_segmentation_rgb")
       error("Wrong sensor selected\n")
    end

    if ~isvarname(varName)
        error("Invalid variable name\n")
    end

    %% Create a python file
    file = fopen(strcat(fileName, '.py'), 'w');

    % Automatically genrates a python file containing the sensor call back
    % bindings
    fprintf(file, 'import numpy as np\n');
    fprintf(file, '\n');
    fprintf(file, '# Converts the tags to human observable colors\n');
    fprintf(file, '_mapDict = { 0:( 0, 0, 0),			1:(70, 70, 70),		2:(190, 153, 153),\n');
    fprintf(file, '             3:(250, 170, 160),		4:(220, 20, 60),	5:(153, 153, 153),\n');
    fprintf(file, '             6:(157, 234, 50),		7:(128, 64, 128),	8:(244, 35, 232),\n');
    fprintf(file, '             9:(107, 142, 35),		10:(0, 0, 142),		11:(102, 102, 156),\n');
    fprintf(file, '             12:(220, 220, 0)}\n');
    fprintf(file, '\n');
    fprintf(file, 'def _setup_mapping():\n');
    fprintf(file, '    global _mapDict\n');
    fprintf(file, '    global _mapping_array\n');
    fprintf(file, '\n');
    fprintf(file, '    tags = np.array(list(_mapDict.keys()))\n');
    fprintf(file, '    colors = np.array(list(_mapDict.values()))\n');
    fprintf(file, '\n');
    fprintf(file, '    _mapping_array = np.zeros((tags.max()+1, 3), dtype="uint8")\n');
    fprintf(file, '    _mapping_array[tags] = colors\n');
    fprintf(file, '\n');

    fprintf(file, 'def bindSensor(sensor):\n');
    fprintf(file, '    _setup_mapping()\n');
    fprintf(file, '    sensor.listen(lambda _image: do_something(_image))\n');
    fprintf(file, '\n');
    fprintf(file, 'def do_something(_image):\n');
    fprintf(file, '    global %s\n', varName);
    fprintf(file, '    global _mapping_array\n');
    fprintf(file, '\n');
    fprintf(file, '    data = np.frombuffer(_image.raw_data, dtype=np.dtype("uint8"))\n');
    fprintf(file, '\n');
    fprintf(file, '    # Get the red channel which has the tags\n');
    fprintf(file, '    data = np.reshape(data, (_image.height, _image.width, 4))\n');
    fprintf(file, '    data = data[:, :, 2]\n');
    fprintf(file, '\n');
    fprintf(file, '    # Map the tags to their respective RGB colors\n');
    fprintf(file, '    data = _mapping_array[data]\n');
    fprintf(file, '\n');
    fprintf(file, '    # Convert the data into MATLAB cast compatible type\n');
    fprintf(file, '    %s = np.ascontiguousarray(data)\n', varName);

    fclose(file);

end

author: darkscyla

repo: MATLAB-Carla-Interface


text:
classdef MparameterThresholder < handle

    properties (Access = public)

    end

    properties (Access = private)
        mV
    end

    properties (Access = private)
       minLengthInUnitCell
       m0V
       m1V
    end

    methods (Access = public)

        function obj = MparameterThresholder(cParams)
            obj.init(cParams);
            obj.m0V = 0;
            obj.m1V = 1;
        end

        function m = thresh(obj,m)
            obj.mV = m;
            obj.thresholdingWhenTooSmallCell();
            obj.thresholdingLimitCases();
            m = obj.mV;
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.minLengthInUnitCell = cParams.minLengthInUnitCell;
        end

        function itIs = isCellTooSmall(obj)
            t = obj.minLengthInUnitCell;
            itIs = t > 1/2;
        end

        function thresholdingWhenTooSmallCell(obj)
            isSmall = obj.isCellTooSmall();
            m = obj.mV;
            m0 = obj.m0V;
            m1 = obj.m1V;
            mT = m;
            mT(isSmall) = m0 + (m1-m0)*heaviside(m(isSmall)-1/2);
            obj.mV = mT;
        end

        function thresholdingLimitCases(obj)
            obj.thresholdingLeftCellFrame();
            obj.thresholdingRightCellFrame();
        end

        function thresholdingLeftCellFrame(obj)
           m = obj.mV;
           t = obj.minLengthInUnitCell;
           mMin = obj.m0V*ones(size(t));
           mMax = mMin + t;
           obj.mV = obj.thresholdCellFrame(m,mMin,mMax);
        end

        function thresholdingRightCellFrame(obj)
           m = obj.mV;
           t = obj.minLengthInUnitCell;
           mMax = obj.m1V*ones(size(t));
           mMin = mMax - t;
           obj.mV = obj.thresholdCellFrame(m,mMin,mMax);
        end

        function m = thresholdCellFrame(obj,m,mMin,mMax)
           isMin = obj.isMin(m,mMin,mMax);
           isMax = obj.isMax(m,mMin,mMax);
           m(isMin) = mMin(isMin);
           m(isMax) = mMax(isMax);
        end

        function itIs = isMin(obj,m,mMin,mMax)
           isFrame        = obj.isInFrameAndCellNotSmall(m,mMin,mMax);
           isCloseToMmin  = obj.isMcloseToMmin(m,mMin,mMax);
           itIs           = isFrame & isCloseToMmin;
        end

        function itIs = isMax(obj,m,mMin,mMax)
           isFrame        = obj.isInFrameAndCellNotSmall(m,mMin,mMax);
           isCloseToMmax  = ~obj.isMcloseToMmin(m,mMin,mMax);
           itIs           = isFrame & isCloseToMmax;
        end

        function itIs = isInFrameAndCellNotSmall(obj,m,mMin,mMax)
           isInFrame      = mMin < m & m < mMax;
           isCellNotSmall = ~obj.isCellTooSmall();
           itIs = isInFrame & isCellNotSmall;
        end

    end

    methods (Access = private, Static)

        function itIs = isMcloseToMmin(m,mMin,mMax)
           itIs = (m - mMin) < (mMax - mMin)/2;
        end

    end

end


text:
function varargout=texthvc(varargin)
%
% For calling details please see v_texthvc.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_texthvc.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: texthvc.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_texthvc(varargin{:});
else
    v_texthvc(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%% OPTI Toolbox Global Nonlinear Programming Examples
%
% This file contains a number of Global NLP problems and demonstrates how
% to solve them using the OPTI Toolbox. You should read and complete
% BasicUsage.m & NonlinearProgramming.m BEFORE running the below examples.
%
%   Copyright (C) 2014 Jonathan Currie (IPL)

% There is also a page on the Wiki which supplements this example:
web('https://www.inverseproblem.co.nz/OPTI/index.php/Probs/GNLP');

%% Determing which Solver to Use
% OPTI Toolbox comes with a number of NLP solvers, thus to determine which
% ones are available on your system you can type:
clc
optiSolver('NLP')

% Note the columns DR and GL. A cross in DR indicates the solver requires
% 1st (and perhaps 2nd) derivatives, while a cross in GL indicates the
% solver can solve Global Optimization problems. For noisy problems unless
% you have exact derivatives, avoid solvers with a cross in DR.

%% Typical Global Optimization Problems
% Global optimization problems result from any of the following circumstances:
%
%       - Objectives containing noise / a stochastic element
%       - Non-convex functions (functions which are not a bowl / hill in 2D)
%       - Objectives that include periodic functions (sin, cos)
%       - Parameter estimation of ODEs solved with adaptive step integrators
%       - Any problem that contains multiple local minima.
%
% An extreme example from Wolfram is shown below:
clc
%Objective
fun = @(x) norm([x(1) - sin(2*x(1) + 3*x(2)) - cos(3*x(1) - 5*x(2));
          x(2) - sin(x(1) - 2*x(2)) + cos(x(1) + 3*x(2))]);
lb = [-4;-4]; ub = [4;4];
x0 = [-4;-4];

%Plot
n = 1e2;
x = linspace(-4,4,n); y = linspace(-4,4,n); Z = zeros(n,n);
for i = 1:n
    for j = 1:n
        Z(j,i) = fun([x(i),y(j)]);
    end
end
surfc(x,y,Z)
colormap summer; shading interp; lighting phong; view(-38,58);
xlabel('x1'); ylabel('x2'); zlabel('obj'); title('Wolfram Global Optimization Problem');

%% Example 1 - Basic Setup
% The main difference when solving a Global NLP is that OPTI treats Global
% and Local NLP problems identically, and therefore you will have to specify
% a Global solver. For this example we will build 4 OPTI objects with 4
% different solvers:
clc
%Build OPTI Problem
prob = optiprob('fun',fun,'bounds',lb,ub);
%Choose Global Solver
opts1 = optiset('solver','nomad');
opts2 = optiset('solver','pswarm');
opts3 = optiset('solver','nlopt','solverOpts',nloptset('algorithm','GN_DIRECT'));
opts4 = optiset('solver','ipopt','warnings','off');

%Pass to OPTI Constructor for Error Checking + Setup
Opt1 = opti(prob,opts1);
Opt2 = opti(prob,opts2);
Opt3 = opti(prob,opts3);
Opt4 = opti(prob,opts4);

% Call solve to solve the problem. Check the plot for a comparison of the
% solution points. Note NOMAD and NLOPT are deterministic with the current
% settings, PSWARM includes random elements, and IPOPT is for comparison of
% a local solution.

[x1,fval1] = solve(Opt1,x0);
[x2,fval2] = solve(Opt2,x0);
[x3,fval3] = solve(Opt3,x0);
[x4,fval4] = solve(Opt4,x0);

view(0,90); hold on;
plot3(x0(1),x0(2),10,'kx','markersize',10);
plot3(x1(1),x1(2),10,'ro'); text(x1(1)+0.1,x1(2)+0.1,10,sprintf('NOMAD: %f',fval1));
plot3(x2(1),x2(2),10,'ro'); text(x2(1)+0.1,x2(2)-0.1,10,sprintf('PSWARM: %f',fval2));
plot3(x3(1),x3(2),10,'ro'); text(x3(1)+0.1,x3(2)+0.2,10,sprintf('NLOPT: %f',fval3));
plot3(x4(1),x4(2),10,'ro'); text(x4(1)+0.1,x4(2)+0.1,10,sprintf('IPOPT: %f',fval4));
hold off;

%% Problem 2 - Quartic
% Includes Linear and Nonlinear Constraints
clc
%Problem
fun = @(x) x(1)^4 - 14*x(1)^2 + 24*x(1) - x(2)^2;

%Linear Constraints
A = [-1 1]; b = 8;
%Nonlinear Constraints
nlcon = @(x) (-x(1)^2) - 2*x(1) + x(2);
nlrhs = -2;
nle = -1;
%Bounds + Starting Guess
lb = [-8;0];
ub = [10;10];
x0 = [0;0];

% Solving a constrained global optimization problem. Note linear
% constraints will be converted to nonlinear ones for this solver
opts = optiset('solver','nomad','solverOpts',nomadset('direction_type','lt 2n'));
Opt = opti('fun',fun,'ineq',A,b,'nlmix',nlcon,nlrhs,nle,'bounds',lb,ub,'options',opts)

% This may take a few seconds...
[x,fval,exitflag,info] = solve(Opt,x0)

%% Example 3 - Solving with PSwarm
% Note x0 is on the local minima side of the saddle
clc
%Problem
fun = @(x) -2*x(1)*x(2);

%Constraints
lb = [-0.5;-0.5];
ub = [1;1];
x0 = [-0.3;-0.3];

%Plot
n = 1e2;
x = linspace(-0.5,1,n); y = linspace(-0.5,1,n); Z = zeros(n,n);
for i = 1:n
    for j = 1:n
        Z(j,i) = fun([x(i),y(j)]);
    end
end
surfc(x,y,Z); hold on; plot3(x0(1),x0(2),fun(x0),'r.','markersize',20); hold off;
colormap winter; shading flat; lighting gouraud; view(18,28);
xlabel('x1'); ylabel('x2'); zlabel('obj'); title('Saddle Point Optimization Problem');

% PSwarm solves bounded and linearly constrained global problems
Opt = opti('fun',fun,'bounds',lb,ub,'options',optiset('solver','pswarm'))

[x,fval,exitflag,info] = solve(Opt,x0)
plot(Opt,3)

%% Example 4 - White Box Quartic
% The following problem is the same as problem 2, however this time we are
% going to solve it using a white box solver (SCIP). The SCIP interface
% will parse the following functions into an algebraic description,
% allowing SCIP to find a global solution to this problem.
clc
%Problem
fun = @(x) x(1)^4 - 14*x(1)^2 + 24*x(1) - x(2)^2;

%Linear Constraints
A = [-1 1]; b = 8;
%Nonlinear Constraints
nlcon = @(x) (-x(1)^2) - 2*x(1) + x(2);
nlrhs = -2;
nle = -1;
%Bounds + Starting Guess
lb = [-8;0];
ub = [10;10];
x0 = [0;0];

% SCIP solves a subset of nonlinear and mixed integer problems, provided
% the problem is deterministics and constains a subset of allowable functions.
Opt = opti('fun',fun,'ineq',A,b,'nlmix',nlcon,nlrhs,nle,'bounds',lb,ub,...
           'options',optiset('solver','scip'))

[x,fval,exitflag,info] = solve(Opt,x0)
plot(Opt)

%% Example 5 - OPTI's Multi-Start Solver
% For simple, low dimensional problems, OPTI's mulit-start is a naive
% algorithm for searching the problem space for a global solution. See the
% following page for more examples:

web('https://www.inverseproblem.co.nz/OPTI/index.php/Advanced/MultiSolve');

clc
Opt = opti(Opt,'solver','ipopt');
[x,fval,exitflag,info] = multisolve(Opt,x0)
multiplot(Opt)

%% Summary
% While Global Optimization solvers may take longer, many real engineering
% problems result in noisy or non-convex objectives and local solvers will
% often struggle to return a result, or fall into the closest local
% optima. OPTI provides a range of competitive blackbox global optimization
% solvers which can return much better results, as well as a white box
% solver for academic users which guarantees a global solution.

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% this is a very naiive and non-optimized cardiac cine GRE sequence

seq=mr.Sequence();              % Create a new sequence object
fov=256e-3; Nx=128; Ny=Nx;      % Define FOV and resolution
alpha=5;                        % flip angle
sliceThickness=5e-3;            % slice
%TE=[7.38 9.84]*1e-3;            % give a vector here to have multiple TEs (e.g. for field mapping)
TE=4.92e-3;
TR=9e-3;                        % only a single value for now

% cardiac features
phases = 8;
hearbeats = 15; % odd numbers of heartbeats / segments work better

% more in-depth parameters
rfSpoilingInc=117;              % RF spoiling increment
rf_duration = 2e-3;
adc_duration = 3.2e-3;
pre_duration = 1e-3;

% set system limits
sys = mr.opts('MaxGrad', 28, 'GradUnit', 'mT/m', ...
    'MaxSlew', 150, 'SlewUnit', 'T/m/s', 'rfRingdownTime', 20e-6, ...
    'rfDeadTime', 100e-6, 'adcDeadTime', 10e-6);

% Create fat-sat pulse
% (in Siemens interpreter from January 2019 duration is limited to 8.192 ms, and although product EPI uses 10.24 ms, 8 ms seems to be sufficient)
% B0=2.89; % 1.5 2.89 3.0
% sat_ppm=-3.45;
% sat_freq=sat_ppm*1e-6*B0*lims.gamma;
% rf_fs = mr.makeGaussPulse(110*pi/180,'system',lims,'Duration',8e-3,...
%     'bandwidth',abs(sat_freq),'freqOffset',sat_freq);
% gz_fs = mr.makeTrapezoid('z',sys,'delay',mr.calcDuration(rf_fs),'Area',1/1e-4); % spoil up to 0.1mm

% define the trigger to play out
trig=mr.makeTrigger('physio1','duration', 2000e-6); % duration after
%trig=mr.makeTriggerPulse('osc0','duration', 4100e-6); % possible channels: 'osc0','osc1','ext1'
trig_out=mr.makeDigitalOutputPulse('ext1','duration', 100e-6,'delay',500e-6); % possible channels: 'osc0','osc1','ext1'

% Create alpha-degree slice selection pulse and gradient
[rf, gz] = mr.makeSincPulse(alpha*pi/180,'Duration',rf_duration,...
    'SliceThickness',sliceThickness,'apodization',0.5,'timeBwProduct',4,'system',sys);

% Define other gradients and ADC events
deltak=1/fov;
gx = mr.makeTrapezoid('x','FlatArea',Nx*deltak,'FlatTime',adc_duration,'system',sys);
adc = mr.makeAdc(Nx,'Duration',gx.flatTime,'Delay',gx.riseTime,'system',sys);
gxPre = mr.makeTrapezoid('x','Area',-gx.area/2,'Duration',pre_duration,'system',sys);
gzReph = mr.makeTrapezoid('z','Area',-gz.area/2,'Duration',pre_duration,'system',sys);

lines_per_segment = round(Ny/hearbeats);
Ns=ceil(Ny/lines_per_segment);
Ny=Ns*lines_per_segment; % it can be that because of the rounding above we measure few more k-space lines...
phaseAreas = ((0:Ny-1)-Ny/2)*deltak;
% now reverse the order in every second segment
phaseAreasSeg=reshape(phaseAreas,lines_per_segment,Ns);
phaseAreasSeg(:,2:2:end)=phaseAreasSeg(end:-1:1,2:2:end);
phaseAreas=phaseAreasSeg(:);

% gradient spoiling
gxSpoil=mr.makeTrapezoid('x','Area',2*Nx*deltak,'system',sys);
gzSpoil=mr.makeTrapezoid('z','Area',4/sliceThickness,'system',sys);

% Calculate timing
delayTE=ceil((TE - mr.calcDuration(gxPre) - gz.fallTime - gz.flatTime/2 ...
    - mr.calcDuration(gx)/2)/seq.gradRasterTime)*seq.gradRasterTime;
delayTR=ceil((TR - mr.calcDuration(gxPre) - mr.calcDuration(gz) ...
    - mr.calcDuration(gx) - delayTE)/seq.gradRasterTime)*seq.gradRasterTime;
assert(all(delayTR>=mr.calcDuration(gxSpoil,gzSpoil)));

fprintf('the sequence will acquire %d lines per segment resulting in a temporal resolution of %g ms per phase\n', lines_per_segment, TR*lines_per_segment*1e3);
fprintf('cardiac acquisition window is: %g ms\n', TR*phases*lines_per_segment*1e3);

rf_phase=0;
rf_inc=0;

% Loop over phase encodes and define sequence blocks
for s=1:Ns
    seq.addBlock(trig); % wait for the cardiac trigger
    for p=1:phases
        for l=1:lines_per_segment
            % restore the line counter
            i=(s-1)*lines_per_segment+l;
            %seq.addBlock(rf_fs,gz_fs); % fat-sat
            rf.phaseOffset=rf_phase/180*pi;
            adc.phaseOffset=rf_phase/180*pi;
            rf_inc=mod(rf_inc+rfSpoilingInc, 360.0);
            rf_phase=mod(rf_phase+rf_inc, 360.0);
            %
            seq.addBlock(rf,gz,trig_out);
            gyPre = mr.makeTrapezoid('y','Area',phaseAreas(i),'Duration',pre_duration,'system',sys);
            seq.addBlock(gxPre,gyPre,gzReph);
            if delayTE>0
                seq.addBlock(mr.makeDelay(delayTE));
            end
            seq.addBlock(gx,adc);
            gyPre.amplitude=-gyPre.amplitude;
            seq.addBlock(mr.makeDelay(delayTR),gxSpoil,gyPre,gzSpoil)
        end
    end
end

%% check whether the timing of the sequence is correct
[ok, error_report]=seq.checkTiming;

if (ok)
    fprintf('Timing check passed successfully\n');
else
    fprintf('Timing check failed! Error listing follows:\n');
    fprintf([error_report{:}]);
    fprintf('\n');
end

%% prepare sequence export
seq.setDefinition('FOV', [fov fov sliceThickness]);
seq.setDefinition('Name', 'cine-gre');

seq.write('cine_gre.seq')       % Write to pulseq file

%seq.install('siemens');
return

%% plot sequence and k-space diagrams

seq.plot('timeRange', [0 5*TR]);

% new single-function call for trajectory calculation
[ktraj_adc, ktraj, t_excitation, t_refocusing, t_adc] = seq.calculateKspace();

% plot k-spaces
time_axis=(1:(size(ktraj,2)))*sys.gradRasterTime;
figure; plot(time_axis, ktraj'); % plot the entire k-space trajectory
hold; plot(t_adc,ktraj_adc(1,:),'.'); % and sampling points on the kx-axis
figure; plot(ktraj(1,:),ktraj(2,:),'b'); % a 2D plot
axis('equal'); % enforce aspect ratio for the correct trajectory display
hold;plot(ktraj_adc(1,:),ktraj_adc(2,:),'r.'); % plot the sampling points

%% very optional slow step, but useful for testing during development e.g. for the real TE, TR or for staying within slewrate limits

rep = seq.testReport;
fprintf([rep{:}]);

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function fixed = p05_fixed_points ( m, fixed_num )

%*****************************************************************************80
%
%% P05_FIXED_POINTS returns the fixed points in problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer FIXED_NUM, the number of fixed points.
%
%    Output, real FIXED(M,FIXED_NUM), the fixed points.
%
  center1 = [  0.0, 0.0 ];
  center2 = [ -0.4, 0.0 ];
  r1 = 1.0;
  r2 = 0.55;

  fixed = [ ...
    center1(1) - r1,  center1(2); ...
    center2(1) - r2,  center2(2); ...
    center2(1) + r2,  center2(2); ...
    center1(1) + r1,  center1(2) ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%PAY Joint forces from payload for SerialLink objects
%
% Calculates the joint loads due to a payload for SerialLink objects,
% It uses the formula Q = J'w, where w is a wrench vector applied at
% the end effector, w = [Fx Fy Fz Mxx Myy Mzz]'. The Jacobian can be
% supplied or computed by RTB
%
% Copyright (C) Bryan Moutrie, 2013-2014
% Licensed under the GNU Lesser General Public License
% see full file for full statement
%
% This file requires file(s) from The Robotics Toolbox for MATLAB (RTB)
% by Peter Corke (www.petercorke.com), see file for statement
%
% Syntax:
%  (1) tauP = pay(w, J)
%  (2) tauP = robot.pay(q, w, f)
%
%  (1) Uses a supplied Jacobian
%  (2) Calculates the Jacobian for joint configuration q in frame f,
%       using either robot.jacob0(q) or robot.jacob0(q) for Jacobian
%
% Outputs:
%  tauP : Generalised joint force/torques
%
% Inputs:
%  robot : SerialLink object with n joints
%  w     : Wrench vector [Fx Fy Fz Mxx Myy Mzz]' in same frame as J
%  J     : Jacobian (supplied): 6-by-n or 6-by-n-m for trajectory
%  q     : Joint row vector, or trajectory matrix of joint row vectors
%  f     : '0' for world frame, 'n' for end-effector frame
%
% See also jacob0, jacobn, paycap, SerialLink.payload

% LICENSE STATEMENT:
%
% This file is part of pHRIWARE.
%
% pHRIWARE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation, either version 3 of
% the License, or (at your option) any later version.
%
% pHRIWARE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public
% License along with pHRIWARE.  If not, see <http://www.gnu.org/licenses/>.
%
% RTB LIBRARY:
%
% Copyright (C) 1993-2014, by Peter I. Corke
% http://www.petercorke.com
% Released under the GNU Lesser General Public license

function tauP = pay(varargin)

if length(varargin) == 3
    w = varargin{2};
    J = varargin{3};
    n = size(J,2);
elseif length(varargin) == 4
    robot = varargin{1};
    q = varargin{2};
    w = varargin{3};
    f = varargin{4};
    n = robot.n;
    J = zeros(6,n,size(q,1));
    if f == '0'
        for i= 1: size(q,1)
            J(:,:,i) = robot.jacob0(q(i,:));
        end
    elseif f == 'n'
        for i= 1: size(q,1)
            J(:,:,i) = robot.jacobn(q(i,:));
        end
    end
end

tauP = -reshape(J(:,:)'*w,n,[])';

end

author: Allopart

repo: rbpf-gmapping


text:
function aid_test ( )

%*****************************************************************************80
%
%% AID_TEST tests R4_AID and R8_AID.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'AID_TEST:\n' );
  fprintf ( 1, '  Test AIRY_AI_PRIME_VALUES, R4_AID, R8_AID\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X     AIRY_AID(X)\n' );
  fprintf ( 1, '                     R4_AID(X)        Diff\n' );
  fprintf ( 1, '                     R8_AID(X)        Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = airy_ai_prime_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_aid ( single ( x ) );
    fx3 = r8_aid ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.6g\n', x, fx1 );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end