text:
function [flips,scorepath,covmats_unflipped] = findflip(data,T,options)
% Finds an optimal permutation of the channels, where goodness is measured
% as the mean lagged partial  cross-correlation across pair of channels and lags.
% In other words, it finds a permutation where the signs of the lagged
% partial correlations agree as much as possible across subjects.
%
% INPUTS
% data          observations, either a struct with X (time series)
%                             or just a matrix containing the time series
%               Alternatively, it can be an (unflipped) array of
%                   autocorrelation matrices (ndim x ndim x no.lags x no. trials),
%                   as computed for example by getCovMats()
% T             length of series. If data is supplied as the
%               autocorrelation matrices, then this is not required and
%               can be specified as []
% options:
%  maxlag        max lag to consider
%  nbatch        no. of channels to evaluate at each iteration (0 for all)
%  noruns        how many random initialisations will be carried out
%  standardise   if 1, standardise the data
%  partial       if 1, base on partial correlation instead of correlation
%  maxcyc        for each initialization, maximum number of cycles of the greedy algorithm
%  verbose       do we get loud?
%
% data_ref and T_ref refer to having a previous data set to
% with which we want to be consistent (see above). These are assumed to be
% already sign-disambiguated

% OUTPUT
% flips         (length(T) X No. channels) binary matrix saying which channels must be
%               flipped for each time series
% scorepath     cell with the score of the winning solutions
% covmats_unflipped  the disambiguated covariance matrices
%
% Author: Diego Vidaurre, University of Oxford.

if nargin < 3, options = struct; end

options = checkoptions_flip(options);
covmats_unflipped = Get_Global_Variables_For_BitFlip_Evaluation(data,T,options);
N = size(covmats_unflipped,4);

ndim = size(covmats_unflipped,1);
score = -Inf;
scorepath = cell(options.noruns,1);

if isinf(options.maxcyc) && (options.nbatch > 0)
    error('If maxcyc is Inf, options.nbatch must be 0')
end

% repetitions of the search
for r = 1:options.noruns

    flipsr = Init_Solution(N,ndim,options,r);
    score_r = EvaluateFlips(flipsr,covmats_unflipped);
    scorepath{r} = score_r;
    if options.verbose
        fprintf('Run %d, Init, score %f \n',r,score_r)
    end

    for cyc = 1:options.maxcyc

        if options.nbatch > 0
            channels = randperm(ndim,options.nbatch);
        else
            channels = randperm(ndim);
        end

        ScoreMatrix = zeros(N,ndim);
        for d = channels
            for j = 1:N
                ScoreMatrix(j,d) = EvaluateFlips(flipsr,covmats_unflipped,j,d);
            end
        end

        [score_r,I] = max(ScoreMatrix(:));
        [j,d] = ind2sub([N ndim],I);

        ds = score_r - max(scorepath{r});
        if ds > options.threshold
            flipsr(j,d) = ~flipsr(j,d);
            if options.verbose
                fprintf('Run %d, Cycle %d, score +%f, flipped (%d,%d) \n',r,cyc,ds,j,d)
            end
            scorepath{r} = [scorepath{r} score_r];
        elseif options.nbatch == 0
            break
        else
            fprintf('Run %d, Cycle %d, score +0 \n',r,cyc)
        end

    end

    if options.verbose
        fprintf('Run %d, Finish, score %f \n',r,scorepath{r}(end))
    end

    if scorepath{r}(end) > score
        score = scorepath{r}(end);
        flips = flipsr;
    end

end

if options.verbose
    fprintf('Final Score=%f\n',score)
end

% Among the equivalent flippings, we keep the one w/ the lowest no. of flips
for j = 1:N
    if mean(flips(j,:))>0.5
        flips(j,:) = 1 - flips(j,:);
    end
end

end

function options = checkoptions_flip(options)
if ~isfield(options,'maxlag'), options.maxlag = 10; end
if ~isfield(options,'noruns'), options.noruns = 1; end
if ~isfield(options,'probinitflip'), options.probinitflip = 0.25; end
if ~isfield(options,'standardise'), options.standardise = 1; end
if ~isfield(options,'partial'), options.partial = 0; end
if ~isfield(options,'verbose'), options.verbose = 1; end
if ~isfield(options,'nbatch'), options.nbatch = 0; end
if ~isfield(options,'maxcyc'), options.maxcyc = 10000; end
if ~isfield(options,'threshold'), options.threshold = 0.00001; end
end

function flipsr = Init_Solution(N,ndim,options,r)
% Prepare the necessary data structures
if ~isfield(options,'Flips')
    if r==1 % first run starts at no flips
        flipsr = zeros(N,ndim);
    else
        flipsr = binornd(1,options.probinitflip,N,ndim); % random init
    end
else
    flipsr = options.Flips;
end
end

function covmats_unflipped = Get_Global_Variables_For_BitFlip_Evaluation(data,T,options)
% Prepare the necessary data structures
if length(size(data))==4 % it is an array of autocorrelation matrices already
    covmats_unflipped = data; clear data
else
    covmats_unflipped = getAllCovMats(data,T,options);
end
end

function s = EvaluateFlips(flipsr,covmats_unflipped,j,d)
% evaluate for a change in subject j and channel d
if nargin>3, flipsr(j,d) = ~flipsr(j,d);  end
signmats = getSignMat(flipsr);
covmats = applySign(covmats_unflipped,signmats);
s = getscore(covmats);
end

function score = getscore(M)
% get the score from the (flipped) matrices of autocorrelation contained in M
% M is a matrix (ndim x ndim x lags x subjects) with autocovariances matrix for all subjects
N = size(M,4); ndim = size(M,1); L = size(M,3);
M = reshape(M,[(ndim^2)*L N]);
C = corr(M);
score = mean(C(triu(true(N),1)));
end

function CovMats = applySign(CovMats,SignMats)
% apply matrices of sign flipping to the autocorrelation matrices
N = size(SignMats,3); nlags = size(CovMats,3);
for j = 1:N
    CovMats(:,:,:,j) = CovMats(:,:,:,j) .* repmat(SignMats(:,:,j),[1 1 nlags]);
end
end

function SignMats = getSignMat(Flips)
% Construct matrices of sign flipping for the autocorrelation matrices
[N,ndim] = size(Flips);
SignMats = zeros(ndim,ndim,N);
for j = 1:N
    flips = ones(1,ndim);
    flips(Flips(j,:)==1) = -1;
    SignMats(:,:,j) = flips' * flips;
end
end

author: OHBA-analysis

repo: HMM-MAR


text:
function gf=filterbankresponse(g,a,L,varargin)
%FILTERBANKRESPONSE  Response of filterbank as function of frequency
%   Usage:  gf=filterbankresponse(g,a,L);
%
%   `gf=filterbankresponse(g,a,L)` computes the total response in frequency
%   of a filterbank specified by *g* and *a* for a signal length of
%   *L*. This corresponds to summing up all channels. The output is a
%   usefull tool to investigate the behaviour of the windows, as peaks
%   indicate that a frequency is overrepresented in the filterbank, while
%   a dip indicates that it is not well represented.
%
%   CAUTION: This function computes a sum of squares of modulus of the
%   frequency responses, which is  also the diagonal of the Fourier
%   transform of the frame operator.
%   Use |filterbankfreqz| for evaluation or plotting of frequency responses
%   of filters.
%
%   `filterbankresponse(g,a,L,'real')` does the same for a filterbank
%   intended for positive-only filterbank.
%
%   `filterbankresponse(g,a,L,fs)` specifies the sampling rate *fs*. This
%   is only used for plotting purposes.
%
%   `gf=filterbankresponse(g,a,L,'individual')` returns responses
%   in frequency of individual filters as columns of a matrix. The total
%   response can be obtained by `gf = sum(gf,2)`.
%
%   `filterbankresponse` takes the following optional parameters:
%
%      'fs',fs
%                 Sampling rate, used only for plotting.
%
%      'complex'
%                 Assume that the filters cover the entire frequency
%                 range. This is the default.
%
%      'real'
%                 Assume that the filters only cover the positive
%                 frequencies (and is intended to work with real-valued
%                 signals only).
%
%      'noplot'
%                 Don't plot the response, just return it.
%
%      'plot'
%                 Plot the response using |plotfftreal| or |plotfft|.
%
%   See also: filterbank, filterbankbounds

definput.flags.ctype={'complex','real'};
definput.flags.plottype={'noplot','plot'};
definput.flags.type={'total','individual'};
definput.keyvals.fs=[];
[flags,kv,fs]=ltfatarghelper({'fs'},definput,varargin);

[g,asan]=filterbankwin(g,a,L,'normal');
M=numel(g);

gf = zeros(L,M);
for m=1:M
    gf(:,m) = comp_filterbankresponse(g(m),asan(m,:),L,flags.do_real);
end

if flags.do_total
    gf = sum(gf,2);
end

if flags.do_plot
    if flags.do_real
        plotfftreal(gf(1:floor(L/2)+1,:),fs,'lin');
    else
        plotfft(gf,fs,'lin');
    end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [x,f,exitflag,output] = minFunc(funObj,x0,options,varargin)
% minFunc(funObj,x0,options,varargin)
%
% Unconstrained optimizer using a line search strategy
%
% Uses an interface very similar to fminunc
%   (it doesn't support all of the optimization toolbox options,
%       but supports many other options).
%
% It computes descent directions using one of ('Method'):
%   - 'sd': Steepest Descent
%       (no previous information used, not recommended)
%   - 'csd': Cyclic Steepest Descent
%       (uses previous step length for a fixed length cycle)
%   - 'bb': Barzilai and Borwein Gradient
%       (uses only previous step)
%   - 'cg': Non-Linear Conjugate Gradient
%       (uses only previous step and a vector beta)
%   - 'scg': Scaled Non-Linear Conjugate Gradient
%       (uses previous step and a vector beta,
%           and Hessian-vector products to initialize line search)
%   - 'pcg': Preconditionined Non-Linear Conjugate Gradient
%       (uses only previous step and a vector beta, preconditioned version)
%   - 'lbfgs': Quasi-Newton with Limited-Memory BFGS Updating
%       (default: uses a predetermined nunber of previous steps to form a
%           low-rank Hessian approximation)
%   - 'newton0': Hessian-Free Newton
%       (numerically computes Hessian-Vector products)
%   - 'pnewton0': Preconditioned Hessian-Free Newton
%       (numerically computes Hessian-Vector products, preconditioned
%       version)
%   - 'qnewton': Quasi-Newton Hessian approximation
%       (uses dense Hessian approximation)
%   - 'mnewton': Newton's method with Hessian calculation after every
%   user-specified number of iterations
%       (needs user-supplied Hessian matrix)
%   - 'newton': Newton's method with Hessian calculation every iteration
%       (needs user-supplied Hessian matrix)
%   - 'tensor': Tensor
%       (needs user-supplied Hessian matrix and Tensor of 3rd partial derivatives)
%
% Several line search strategies are available for finding a step length satisfying
%   the termination criteria ('LS'):
%   - 0: Backtrack w/ Step Size Halving
%   - 1: Backtrack w/ Quadratic/Cubic Interpolation from new function values
%   - 2: Backtrack w/ Cubic Interpolation from new function + gradient
%   values (default for 'bb' and 'sd')
%   - 3: Bracketing w/ Step Size Doubling and Bisection
%   - 4: Bracketing w/ Cubic Interpolation/Extrapolation with function +
%   gradient values (default for all except 'bb' and 'sd')
%   - 5: Bracketing w/ Mixed Quadratic/Cubic Interpolation/Extrapolation
%   - 6: Use Matlab Optimization Toolbox's line search
%           (requires Matlab's linesearch.m to be added to the path)
%
%   Above, the first three find a point satisfying the Armijo conditions,
%   while the last four search for find a point satisfying the Wolfe
%   conditions.  If the objective function overflows, it is recommended
%   to use one of the first 3.
%   The first three can be used to perform a non-monotone
%   linesearch by changing the option 'Fref'.
%
% Several strategies for choosing the initial step size are avaiable ('LS_init'):
%   - 0: Always try an initial step length of 1 (default for all except 'cg' and 'sd')
%       (t = 1)
%   - 1: Use a step similar to the previous step (default for 'cg' and 'sd')
%       (t = t_old*min(2,g'd/g_old'd_old))
%   - 2: Quadratic Initialization using previous function value and new
%   function value/gradient (use this if steps tend to be very long)
%       (t = min(1,2*(f-f_old)/g))
%   - 3: The minimum between 1 and twice the previous step length
%       (t = min(1,2*t)
%   - 4: The scaled conjugate gradient step length (may accelerate
%   conjugate gradient methods, but requires a Hessian-vector product)
%       (t = g'd/d'Hd)
%
% Inputs:
%   funObj is a function handle
%   x0 is a starting vector;
%   options is a struct containing parameters
%  (defaults are used for non-existent or blank fields)
%   all other arguments are passed to funObj
%
% Outputs:
%   x is the minimum value found
%   f is the function value at the minimum found
%   exitflag returns an exit condition
%   output returns a structure with other information
%
% Supported Input Options
%   Display - Level of display [ off | final | (iter) | full | excessive ]
%   MaxFunEvals - Maximum number of function evaluations allowed (1000)
%   MaxIter - Maximum number of iterations allowed (500)
%   TolFun - Termination tolerance on the first-order optimality (1e-5)
%   TolX - Termination tolerance on progress in terms of function/parameter changes (1e-9)
%   Method - [ sd | csd | bb | cg | scg | pcg | {lbfgs} | newton0 | pnewton0 |
%       qnewton | mnewton | newton | tensor ]
%   c1 - Sufficient Decrease for Armijo condition (1e-4)
%   c2 - Curvature Decrease for Wolfe conditions (.2 for cg methods, .9 otherwise)
%   LS_init - Line Search Initialization -see above (2 for cg/sd, 4 for scg, 0 otherwise)
%   LS - Line Search type -see above (2 for bb, 4 otherwise)
%   Fref - Setting this to a positive integer greater than 1
%       will use non-monotone Armijo objective in the line search.
%       (20 for bb, 10 for csd, 1 for all others)
%   numDiff - compute derivative numerically
%       (default: 0) (this option has a different effect for 'newton', see below)
%   useComplex - if 1, use complex differentials when computing numerical derivatives
%       to get very accurate values (default: 0)
%   DerivativeCheck - if 'on', computes derivatives numerically at initial
%       point and compares to user-supplied derivative (default: 'off')
%   outputFcn - function to run after each iteration (default: []).  It
%       should have the following interface:
%       outputFcn(x,infoStruct,state,varargin{:})
%   useMex - where applicable, use mex files to speed things up (default: 1)
%
% Method-specific input options:
%   newton:
%       HessianModify - type of Hessian modification for direct solvers to
%       use if the Hessian is not positive definite (default: 0)
%           0: Minimum Euclidean norm s.t. eigenvalues sufficiently large
%           (requires eigenvalues on iterations where matrix is not pd)
%           1: Start with (1/2)*||A||_F and increment until Cholesky succeeds
%           (an approximation to method 0, does not require eigenvalues)
%           2: Modified LDL factorization
%           (only 1 generalized Cholesky factorization done and no eigenvalues required)
%           3: Modified Spectral Decomposition
%           (requires eigenvalues)
%           4: Modified Symmetric Indefinite Factorization
%           5: Uses the eigenvector of the smallest eigenvalue as negative
%           curvature direction
%       cgSolve - use conjugate gradient instead of direct solver (default: 0)
%           0: Direct Solver
%           1: Conjugate Gradient
%           2: Conjugate Gradient with Diagonal Preconditioner
%           3: Conjugate Gradient with LBFGS Preconditioner
%           x: Conjugate Graident with Symmetric Successive Over Relaxation
%           Preconditioner with parameter x
%               (where x is a real number in the range [0,2])
%           x: Conjugate Gradient with Incomplete Cholesky Preconditioner
%           with drop tolerance -x
%               (where x is a real negative number)
%       numDiff - compute Hessian numerically
%                 (default: 0, done with complex differentials if useComplex = 1)
%       LS_saveHessiancomp - when on, only computes the Hessian at the
%       first and last iteration of the line search (default: 1)
%   mnewton:
%       HessianIter - number of iterations to use same Hessian (default: 5)
%   qnewton:
%       initialHessType - scale initial Hessian approximation (default: 1)
%       qnUpdate - type of quasi-Newton update (default: 3):
%           0: BFGS
%           1: SR1 (when it is positive-definite, otherwise BFGS)
%           2: Hoshino
%           3: Self-Scaling BFGS
%           4: Oren's Self-Scaling Variable Metric method
%           5: McCormick-Huang asymmetric update
%       Damped - use damped BFGS update (default: 1)
%   newton0/pnewton0:
%       HvFunc - user-supplied function that returns Hessian-vector products
%           (by default, these are computed numerically using autoHv)
%           HvFunc should have the following interface: HvFunc(v,x,varargin{:})
%       useComplex - use a complex perturbation to get high accuracy
%           Hessian-vector products (default: 0)
%           (the increased accuracy can make the method much more efficient,
%               but gradient code must properly support complex inputs)
%       useNegCurv - a negative curvature direction is used as the descent
%           direction if one is encountered during the cg iterations
%           (default: 1)
%       precFunc (for pnewton0 only) - user-supplied preconditioner
%           (by default, an L-BFGS preconditioner is used)
%           precFunc should have the following interfact:
%           precFunc(v,x,varargin{:})
%   lbfgs:
%       Corr - number of corrections to store in memory (default: 100)
%           (higher numbers converge faster but use more memory)
%       Damped - use damped update (default: 0)
%   pcg:
%       cgUpdate - type of update (default: 2)
%   cg/scg/pcg:
%       cgUpdate - type of update (default for cg/scg: 2, default for pcg: 1)
%           0: Fletcher Reeves
%           1: Polak-Ribiere
%           2: Hestenes-Stiefel (not supported for pcg)
%           3: Gilbert-Nocedal
%       HvFunc (for scg only)- user-supplied function that returns Hessian-vector
%           products
%           (by default, these are computed numerically using autoHv)
%           HvFunc should have the following interface:
%           HvFunc(v,x,varargin{:})
%       precFunc (for pcg only) - user-supplied preconditioner
%           (by default, an L-BFGS preconditioner is used)
%           precFunc should have the following interfact:
%           precFunc(v,x,varargin{:})
%   bb:
%       bbType - type of bb step (default: 1)
%           0: min_alpha ||delta_x - alpha delta_g||_2
%           1: min_alpha ||alpha delta_x - delta_g||_2
%           2: Conic BB
%           3: Gradient method with retards
%   csd:
%       cycle - length of cycle (default: 3)
%
% Supported Output Options
%   iterations - number of iterations taken
%   funcCount - number of function evaluations
%   algorithm - algorithm used
%   firstorderopt - first-order optimality
%   message - exit message
%   trace.funccount - function evaluations after each iteration
%   trace.fval - function value after each iteration
%
% Author: Mark Schmidt (2006)
% Web: http://www.cs.ubc.ca/~schmidtm
%
% Sources (in order of how much the source material contributes):
%   J. Nocedal and S.J. Wright.  1999.  "Numerical Optimization".  Springer Verlag.
%   R. Fletcher.  1987.  "Practical Methods of Optimization".  Wiley.
%   J. Demmel.  1997.  "Applied Linear Algebra.  SIAM.
%   R. Barret, M. Berry, T. Chan, J. Demmel, J. Dongarra, V. Eijkhout, R.
%   Pozo, C. Romine, and H. Van der Vost.  1994.  "Templates for the Solution of
%   Linear Systems: Building Blocks for Iterative Methods".  SIAM.
%   J. More and D. Thuente.  "Line search algorithms with guaranteed
%   sufficient decrease".  ACM Trans. Math. Softw. vol 20, 286-307, 1994.
%   M. Raydan.  "The Barzilai and Borwein gradient method for the large
%   scale unconstrained minimization problem".  SIAM J. Optim., 7, 26-33,
%   (1997).
%   "Mathematical Optimization".  The Computational Science Education
%   Project.  1995.
%   C. Kelley.  1999.  "Iterative Methods for Optimization".  Frontiers in
%   Applied Mathematics.  SIAM.

if nargin < 3
    options = [];
end

% Get Parameters
[verbose,verboseI,debug,doPlot,maxFunEvals,maxIter,tolFun,tolX,method,...
    corrections,c1,c2,LS_init,LS,cgSolve,qnUpdate,cgUpdate,initialHessType,...
    HessianModify,Fref,useComplex,numDiff,LS_saveHessianComp,...
    DerivativeCheck,Damped,HvFunc,bbType,cycle,...
    HessianIter,outputFcn,useMex,useNegCurv,precFunc] = ...
    minFunc_processInputOptions(options);

% Constants
SD = 0;
CSD = 1;
BB = 2;
CG = 3;
PCG = 4;
LBFGS = 5;
QNEWTON = 6;
NEWTON0 = 7;
NEWTON = 8;
TENSOR = 9;

% Initialize
p = length(x0);
d = zeros(p,1);
x = x0;
t = 1;

% If necessary, form numerical differentiation functions
funEvalMultiplier = 1;
if numDiff && method ~= TENSOR
    varargin(3:end+2) = varargin(1:end);
    varargin{1} = useComplex;
    varargin{2} = funObj;
    if method ~= NEWTON
        if debug
            if useComplex
                fprintf('Using complex differentials for gradient computation\n');
            else
                fprintf('Using finite differences for gradient computation\n');
            end
        end
        funObj = @autoGrad;
    else
        if debug
            if useComplex
                fprintf('Using complex differentials for gradient computation\n');
            else
                fprintf('Using finite differences for gradient computation\n');
            end
        end
        funObj = @autoHess;
    end

    if method == NEWTON0 && useComplex == 1
        if debug
            fprintf('Turning off the use of complex differentials\n');
        end
        useComplex = 0;
    end

    if useComplex
        funEvalMultiplier = p;
    else
        funEvalMultiplier = p+1;
    end
end

% Evaluate Initial Point
if method < NEWTON
    [f,g] = funObj(x,varargin{:});
else
    [f,g,H] = funObj(x,varargin{:});
    computeHessian = 1;
end
funEvals = 1;

if strcmp(DerivativeCheck,'on')
    if numDiff
        fprintf('Can not do derivative checking when numDiff is 1\n');
    end
    % Check provided gradient/hessian function using numerical derivatives
    fprintf('Checking Gradient:\n');
    [f2,g2] = autoGrad(x,useComplex,funObj,varargin{:});

    fprintf('Max difference between user and numerical gradient: %f\n',max(abs(g-g2)));
    if max(abs(g-g2)) > 1e-4
        fprintf('User NumDif:\n');
        [g g2]
        diff = abs(g-g2)
        pause;
    end

    if method >= NEWTON
        fprintf('Check Hessian:\n');
        [f2,g2,H2] = autoHess(x,useComplex,funObj,varargin{:});

        fprintf('Max difference between user and numerical hessian: %f\n',max(abs(H(:)-H2(:))));
        if max(abs(H(:)-H2(:))) > 1e-4
            H
            H2
            diff = abs(H-H2)
            pause;
        end
    end
end

% Output Log
if verboseI
    fprintf('%10s %10s %15s %15s %15s\n','Iteration','FunEvals','Step Length','Function Val','Opt Cond');
end

% Output Function
if ~isempty(outputFcn)
    callOutput(outputFcn,x,'init',0,funEvals,f,[],[],g,[],sum(abs(g)),varargin{:});
end

% Initialize Trace
trace.fval = f;
trace.funcCount = funEvals;

% Check optimality of initial point
if sum(abs(g)) <= tolFun
    exitflag=1;
    msg = 'Optimality Condition below TolFun';
    if verbose
        fprintf('%s\n',msg);
    end
    if nargout > 3
        output = struct('iterations',0,'funcCount',1,...
            'algorithm',method,'firstorderopt',sum(abs(g)),'message',msg,'trace',trace);
    end
    return;
end

% Perform up to a maximum of 'maxIter' descent steps:
for i = 1:maxIter

    % ****************** COMPUTE DESCENT DIRECTION *****************

    switch method
        case SD % Steepest Descent
            d = -g;

        case CSD % Cyclic Steepest Descent

            if mod(i,cycle) == 1 % Use Steepest Descent
                alpha = 1;
                LS_init = 2;
                LS = 4; % Precise Line Search
            elseif mod(i,cycle) == mod(1+1,cycle) % Use Previous Step
                alpha = t;
                LS_init = 0;
                LS = 2; % Non-monotonic line search
            end
            d = -alpha*g;

        case BB % Steepest Descent with Barzilai and Borwein Step Length

            if i == 1
                d = -g;
            else
                y = g-g_old;
                s = t*d;
                if bbType == 0
                    yy = y'*y;
                    alpha = (s'*y)/(yy);
                    if alpha <= 1e-10 || alpha > 1e10
                        alpha = 1;
                    end
                elseif bbType == 1
                    sy = s'*y;
                    alpha = (s'*s)/sy;
                    if alpha <= 1e-10 || alpha > 1e10
                        alpha = 1;
                    end
                elseif bbType == 2 % Conic Interpolation ('Modified BB')
                    sy = s'*y;
                    ss = s'*s;
                    alpha = ss/sy;
                    if alpha <= 1e-10 || alpha > 1e10
                        alpha = 1;
                    end
                    alphaConic = ss/(6*(myF_old - f) + 4*g'*s + 2*g_old'*s);
                    if alphaConic > .001*alpha && alphaConic < 1000*alpha
                        alpha = alphaConic;
                    end
                elseif bbType == 3 % Gradient Method with retards (bb type 1, random selection of previous step)
                    sy = s'*y;
                    alpha = (s'*s)/sy;
                    if alpha <= 1e-10 || alpha > 1e10
                        alpha = 1;
                    end
                    v(1+mod(i-2,5)) = alpha;
                    alpha = v(ceil(rand*length(v)));
                end
                d = -alpha*g;
            end
            g_old = g;
            myF_old = f;

        case CG % Non-Linear Conjugate Gradient

            if i == 1
                d = -g; % Initially use steepest descent direction
            else
                gtgo = g'*g_old;
                gotgo = g_old'*g_old;

                if cgUpdate == 0
                    % Fletcher-Reeves
                    beta = (g'*g)/(gotgo);
                elseif cgUpdate == 1
                    % Polak-Ribiere
                    beta = (g'*(g-g_old)) /(gotgo);
                elseif cgUpdate == 2
                    % Hestenes-Stiefel
                    beta = (g'*(g-g_old))/((g-g_old)'*d);
                else
                    % Gilbert-Nocedal
                    beta_FR = (g'*(g-g_old)) /(gotgo);
                    beta_PR = (g'*g-gtgo)/(gotgo);
                    beta = max(-beta_FR,min(beta_PR,beta_FR));
                end

                d = -g + beta*d;

                % Restart if not a direction of sufficient descent
                if g'*d > -tolX
                    if debug
                        fprintf('Restarting CG\n');
                    end
                    beta = 0;
                    d = -g;
                end

                % Old restart rule:
                %if beta < 0 || abs(gtgo)/(gotgo) >= 0.1 || g'*d >= 0

            end
            g_old = g;

        case PCG % Preconditioned Non-Linear Conjugate Gradient

            % Apply preconditioner to negative gradient
            if isempty(precFunc)
                % Use L-BFGS Preconditioner
                if i == 1
                    old_dirs = zeros(length(g),0);
                    old_stps = zeros(length(g),0);
                    Hdiag = 1;
                    s = -g;
                else
                    [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);

                    if useMex
                        s = lbfgsC(-g,old_dirs,old_stps,Hdiag);
                    else
                        s = lbfgs(-g,old_dirs,old_stps,Hdiag);
                    end
                end
            else % User-supplied preconditioner
                s = precFunc(-g,x,varargin{:});
            end

            if i == 1
                d = s;
            else

                if cgUpdate == 0
                    % Preconditioned Fletcher-Reeves
                    beta = (g'*s)/(g_old'*s_old);
                elseif cgUpdate < 3
                    % Preconditioned Polak-Ribiere
                    beta = (g'*(s-s_old))/(g_old'*s_old);
                else
                    % Preconditioned Gilbert-Nocedal
                    beta_FR = (g'*s)/(g_old'*s_old);
                    beta_PR = (g'*(s-s_old))/(g_old'*s_old);
                    beta = max(-beta_FR,min(beta_PR,beta_FR));
                end
                d = s + beta*d;

                if g'*d > -tolX
                    if debug
                        fprintf('Restarting CG\n');
                    end
                    beta = 0;
                    d = s;
                end

            end
            g_old = g;
            s_old = s;
        case LBFGS % L-BFGS

            % Update the direction and step sizes

            if i == 1
                d = -g; % Initially use steepest descent direction
                old_dirs = zeros(length(g),0);
                old_stps = zeros(length(d),0);
                Hdiag = 1;
            else
                if Damped
                    [old_dirs,old_stps,Hdiag] = dampedUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);
                else
                    [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);
                end

                if useMex
                    d = lbfgsC(-g,old_dirs,old_stps,Hdiag);
                else
                    d = lbfgs(-g,old_dirs,old_stps,Hdiag);
                end
            end
            g_old = g;

        case QNEWTON % Use quasi-Newton Hessian approximation

            if i == 1
                d = -g;
            else
                % Compute difference vectors
                y = g-g_old;
                s = t*d;

                if i == 2
                    % Make initial Hessian approximation
                    if initialHessType == 0
                        % Identity
                        if qnUpdate <= 1
                            R = eye(length(g));
                        else
                            H = eye(length(g));
                        end
                    else
                        % Scaled Identity
                        if debug
                            fprintf('Scaling Initial Hessian Approximation\n');
                        end
                        if qnUpdate <= 1
                            % Use Cholesky of Hessian approximation
                            R = sqrt((y'*y)/(y'*s))*eye(length(g));
                        else
                            % Use Inverse of Hessian approximation
                            H = eye(length(g))*(y'*s)/(y'*y);
                        end
                    end
                end

                if qnUpdate == 0 % Use BFGS updates
                    Bs = R'*(R*s);
                    if Damped
                        eta = .02;
                        if y'*s < eta*s'*Bs
                            if debug
                                fprintf('Damped Update\n');
                            end
                            theta = min(max(0,((1-eta)*s'*Bs)/(s'*Bs - y'*s)),1);
                            y = theta*y + (1-theta)*Bs;
                        end
                        R = cholupdate(cholupdate(R,y/sqrt(y'*s)),Bs/sqrt(s'*Bs),'-');
                    else
                        if y'*s > 1e-10
                            R = cholupdate(cholupdate(R,y/sqrt(y'*s)),Bs/sqrt(s'*Bs),'-');
                        else
                            if debug
                                fprintf('Skipping Update\n');
                            end
                        end
                    end
                elseif qnUpdate == 1 % Perform SR1 Update if it maintains positive-definiteness

                    Bs = R'*(R*s);
                    ymBs = y-Bs;
                    if abs(s'*ymBs) >= norm(s)*norm(ymBs)*1e-8 && (s-((R\(R'\y))))'*y > 1e-10
                        R = cholupdate(R,-ymBs/sqrt(ymBs'*s),'-');
                    else
                        if debug
                            fprintf('SR1 not positive-definite, doing BFGS Update\n');
                        end
                        if Damped
                            eta = .02;
                            if y'*s < eta*s'*Bs
                                if debug
                                    fprintf('Damped Update\n');
                                end
                                theta = min(max(0,((1-eta)*s'*Bs)/(s'*Bs - y'*s)),1);
                                y = theta*y + (1-theta)*Bs;
                            end
                            R = cholupdate(cholupdate(R,y/sqrt(y'*s)),Bs/sqrt(s'*Bs),'-');
                        else
                            if y'*s > 1e-10
                                R = cholupdate(cholupdate(R,y/sqrt(y'*s)),Bs/sqrt(s'*Bs),'-');
                            else
                                if debug
                                    fprintf('Skipping Update\n');
                                end
                            end
                        end
                    end
                elseif qnUpdate == 2 % Use Hoshino update
                    v = sqrt(y'*H*y)*(s/(s'*y) - (H*y)/(y'*H*y));
                    phi = 1/(1 + (y'*H*y)/(s'*y));
                    H = H + (s*s')/(s'*y) - (H*y*y'*H)/(y'*H*y) + phi*v*v';

                elseif qnUpdate == 3 % Self-Scaling BFGS update
                    ys = y'*s;
                    Hy = H*y;
                    yHy = y'*Hy;
                    gamma = ys/yHy;
                    v = sqrt(yHy)*(s/ys - Hy/yHy);
                    H = gamma*(H - Hy*Hy'/yHy + v*v') + (s*s')/ys;
                elseif qnUpdate == 4 % Oren's Self-Scaling Variable Metric update

                    % Oren's method
                    if (s'*y)/(y'*H*y) > 1
                        phi = 1; % BFGS
                        omega = 0;
                    elseif (s'*(H\s))/(s'*y) < 1
                        phi = 0; % DFP
                        omega = 1;
                    else
                        phi = (s'*y)*(y'*H*y-s'*y)/((s'*(H\s))*(y'*H*y)-(s'*y)^2);
                        omega = phi;
                    end

                    gamma = (1-omega)*(s'*y)/(y'*H*y) + omega*(s'*(H\s))/(s'*y);
                    v = sqrt(y'*H*y)*(s/(s'*y) - (H*y)/(y'*H*y));
                    H = gamma*(H - (H*y*y'*H)/(y'*H*y) + phi*v*v') + (s*s')/(s'*y);

                elseif qnUpdate == 5 % McCormick-Huang asymmetric update
                    theta = 1;
                    phi = 0;
                    psi = 1;
                    omega = 0;
                    t1 = s*(theta*s + phi*H'*y)';
                    t2 = (theta*s + phi*H'*y)'*y;
                    t3 = H*y*(psi*s + omega*H'*y)';
                    t4 = (psi*s + omega*H'*y)'*y;
                    H = H + t1/t2 - t3/t4;
                end

                if qnUpdate <= 1
                    d = -R\(R'\g);
                else
                    d = -H*g;
                end

            end
            g_old = g;

        case NEWTON0 % Hessian-Free Newton

            cgMaxIter = min(p,maxFunEvals-funEvals);
            cgForce = min(0.5,sqrt(norm(g)))*norm(g);

            % Set-up preconditioner
            precondFunc = [];
            precondArgs = [];
            if cgSolve == 1
                if isempty(precFunc) % Apply L-BFGS preconditioner
                    if i == 1
                        old_dirs = zeros(length(g),0);
                        old_stps = zeros(length(g),0);
                        Hdiag = 1;
                    else
                        [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);
                        if useMex
                            precondFunc = @lbfgsC;
                        else
                            precondFunc = @lbfgs;
                        end
                        precondArgs = {old_dirs,old_stps,Hdiag};
                    end
                    g_old = g;
                else
                    % Apply user-defined preconditioner
                    precondFunc = precFunc;
                    precondArgs = {x,varargin{:}};
                end
            end

            % Solve Newton system using cg and hessian-vector products
            if isempty(HvFunc)
                % No user-supplied Hessian-vector function,
                % use automatic differentiation
                HvFun = @autoHv;
                HvArgs = {x,g,useComplex,funObj,varargin{:}};
            else
                % Use user-supplid Hessian-vector function
                HvFun = HvFunc;
                HvArgs = {x,varargin{:}};
            end

            if useNegCurv
                [d,cgIter,cgRes,negCurv] = conjGrad([],-g,cgForce,cgMaxIter,debug,precondFunc,precondArgs,HvFun,HvArgs);
            else
                [d,cgIter,cgRes] = conjGrad([],-g,cgForce,cgMaxIter,debug,precondFunc,precondArgs,HvFun,HvArgs);
            end

            funEvals = funEvals+cgIter;
            if debug
                fprintf('newtonCG stopped on iteration %d w/ residual %.5e\n',cgIter,cgRes);

            end

            if useNegCurv
                if ~isempty(negCurv)
                    %if debug
                    fprintf('Using negative curvature direction\n');
                    %end
                    d = negCurv/norm(negCurv);
                    d = d/sum(abs(g));
                end
            end

        case NEWTON % Newton search direction

            if cgSolve == 0
                if HessianModify == 0
                    % Attempt to perform a Cholesky factorization of the Hessian
                    [R,posDef] = chol(H);

                    % If the Cholesky factorization was successful, then the Hessian is
                    % positive definite, solve the system
                    if posDef == 0
                        d = -R\(R'\g);

                    else
                        % otherwise, adjust the Hessian to be positive definite based on the
                        % minimum eigenvalue, and solve with QR
                        % (expensive, we don't want to do this very much)
                        if debug
                            fprintf('Adjusting Hessian\n');
                        end
                        H = H + eye(length(g)) * max(0,1e-12 - min(real(eig(H))));
                        d = -H\g;
                    end
                elseif HessianModify == 1
                    % Modified Incomplete Cholesky
                    R = mcholinc(H,debug);
                    d = -R\(R'\g);
                elseif HessianModify == 2
                    % Modified Generalized Cholesky
                    if useMex
                        [L D perm] = mcholC(H);
                    else
                        [L D perm] = mchol(H);
                    end
                    d(perm) = -L' \ ((D.^-1).*(L \ g(perm)));

                elseif HessianModify == 3
                    % Modified Spectral Decomposition
                    [V,D] = eig((H+H')/2);
                    D = diag(D);
                    D = max(abs(D),max(max(abs(D)),1)*1e-12);
                    d = -V*((V'*g)./D);
                elseif HessianModify == 4
                    % Modified Symmetric Indefinite Factorization
                    [L,D,perm] = ldl(H,'vector');
                    [blockPos junk] = find(triu(D,1));
                    for diagInd = setdiff(setdiff(1:p,blockPos),blockPos+1)
                        if D(diagInd,diagInd) < 1e-12
                            D(diagInd,diagInd) = 1e-12;
                        end
                    end
                    for blockInd = blockPos'
                        block = D(blockInd:blockInd+1,blockInd:blockInd+1);
                        block_a = block(1);
                        block_b = block(2);
                        block_d = block(4);
                        lambda = (block_a+block_d)/2 - sqrt(4*block_b^2 + (block_a - block_d)^2)/2;
                        D(blockInd:blockInd+1,blockInd:blockInd+1) = block+eye(2)*(lambda+1e-12);
                    end
                    d(perm) = -L' \ (D \ (L \ g(perm)));
                else
                    % Take Newton step if Hessian is pd,
                    % otherwise take a step with negative curvature
                    [R,posDef] = chol(H);
                    if posDef == 0
                        d = -R\(R'\g);
                    else
                        if debug
                            fprintf('Taking Direction of Negative Curvature\n');
                        end
                        [V,D] = eig(H);
                        u = V(:,1);
                        d = -sign(u'*g)*u;
                    end
                end

            else
                % Solve with Conjugate Gradient
                cgMaxIter = p;
                cgForce = min(0.5,sqrt(norm(g)))*norm(g);

                % Select Preconditioner
                if cgSolve == 1
                    % No preconditioner
                    precondFunc = [];
                    precondArgs = [];
                elseif cgSolve == 2
                    % Diagonal preconditioner
                    precDiag = diag(H);
                    precDiag(precDiag < 1e-12) = 1e-12 - min(precDiag);
                    precondFunc = @precondDiag;
                    precondArgs = {precDiag.^-1};
                elseif cgSolve == 3
                    % L-BFGS preconditioner
                    if i == 1
                        old_dirs = zeros(length(g),0);
                        old_stps = zeros(length(g),0);
                        Hdiag = 1;
                    else
                        [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);
                    end
                    g_old = g;
                    if useMex
                        precondFunc = @lbfgsC;
                    else
                        precondFunc = @lbfgs;
                    end
                    precondArgs = {old_dirs,old_stps,Hdiag};
                elseif cgSolve > 0
                    % Symmetric Successive Overelaxation Preconditioner
                    omega = cgSolve;
                    D = diag(H);
                    D(D < 1e-12) = 1e-12 - min(D);
                    precDiag = (omega/(2-omega))*D.^-1;
                    precTriu = diag(D/omega) + triu(H,1);
                    precondFunc = @precondTriuDiag;
                    precondArgs = {precTriu,precDiag.^-1};
                else
                    % Incomplete Cholesky Preconditioner
                    opts.droptol = -cgSolve;
                    opts.rdiag = 1;
                    R = cholinc(sparse(H),opts);
                    if min(diag(R)) < 1e-12
                        R = cholinc(sparse(H + eye*(1e-12 - min(diag(R)))),opts);
                    end
                    precondFunc = @precondTriu;
                    precondArgs = {R};
                end

                % Run cg with the appropriate preconditioner
                if isempty(HvFunc)
                    % No user-supplied Hessian-vector function
                    [d,cgIter,cgRes] = conjGrad(H,-g,cgForce,cgMaxIter,debug,precondFunc,precondArgs);
                else
                    % Use user-supplied Hessian-vector function
                    [d,cgIter,cgRes] = conjGrad(H,-g,cgForce,cgMaxIter,debug,precondFunc,precondArgs,HvFunc,{x,varargin{:}});
                end
                if debug
                    fprintf('CG stopped after %d iterations w/ residual %.5e\n',cgIter,cgRes);
                    %funEvals = funEvals + cgIter;
                end
            end

        case TENSOR % Tensor Method

            if numDiff
                % Compute 3rd-order Tensor Numerically
                [junk1 junk2 junk3 T] = autoTensor(x,useComplex,funObj,varargin{:});
            else
                % Use user-supplied 3rd-derivative Tensor
                [junk1 junk2 junk3 T] = funObj(x,varargin{:});
            end
            options_sub.Method = 'newton';
            options_sub.Display = 'none';
            options_sub.TolX = tolX;
            options_sub.TolFun = tolFun;
            d = minFunc(@taylorModel,zeros(p,1),options_sub,f,g,H,T);

            if any(abs(d) > 1e5) || all(abs(d) < 1e-5) || g'*d > -tolX
                if debug
                    fprintf('Using 2nd-Order Step\n');
                end
                [V,D] = eig((H+H')/2);
                D = diag(D);
                D = max(abs(D),max(max(abs(D)),1)*1e-12);
                d = -V*((V'*g)./D);
            else
                if debug
                    fprintf('Using 3rd-Order Step\n');
                end
            end
    end

    if ~isLegal(d)
        fprintf('Step direction is illegal!\n');
        pause;
        return
    end

    % ****************** COMPUTE STEP LENGTH ************************

    % Directional Derivative
    gtd = g'*d;

    % Check that progress can be made along direction
    if gtd > -tolX
        exitflag=2;
        msg = 'Directional Derivative below TolX';
        break;
    end

    % Select Initial Guess
    if i == 1
        if method < NEWTON0
            t = min(1,1/sum(abs(g)));
        else
            t = 1;
        end
    else
        if LS_init == 0
            % Newton step
            t = 1;
        elseif LS_init == 1
            % Close to previous step length
            t = t*min(2,(gtd_old)/(gtd));
        elseif LS_init == 2
            % Quadratic Initialization based on {f,g} and previous f
            t = min(1,2*(f-f_old)/(gtd));
        elseif LS_init == 3
            % Double previous step length
            t = min(1,t*2);
        elseif LS_init == 4
            % Scaled step length if possible
            if isempty(HvFunc)
                % No user-supplied Hessian-vector function,
                % use automatic differentiation
                dHd = d'*autoHv(d,x,g,0,funObj,varargin{:});
            else
                % Use user-supplid Hessian-vector function
                dHd = d'*HvFunc(d,x,varargin{:});
            end

            funEvals = funEvals + 1;
            if dHd > 0
                t = -gtd/(dHd);
            else
                t = min(1,2*(f-f_old)/(gtd));
            end
        end

        if t <= 0
            t = 1;
        end
    end
    f_old = f;
    gtd_old = gtd;

    % Compute reference fr if using non-monotone objective
    if Fref == 1
        fr = f;
    else
        if i == 1
            old_fvals = repmat(-inf,[Fref 1]);
        end

        if i <= Fref
            old_fvals(i) = f;
        else
            old_fvals = [old_fvals(2:end);f];
        end
        fr = max(old_fvals);
    end

    computeHessian = 0;
    if method >= NEWTON
        if HessianIter == 1
            computeHessian = 1;
        elseif i > 1 && mod(i-1,HessianIter) == 0
            computeHessian = 1;
        end
    end

    % Line Search
    f_old = f;
    if LS < 3 % Use Armijo Bactracking
        % Perform Backtracking line search
        if computeHessian
            [t,x,f,g,LSfunEvals,H] = ArmijoBacktrack(x,t,d,f,fr,g,gtd,c1,LS,tolX,debug,doPlot,LS_saveHessianComp,funObj,varargin{:});
        else
            [t,x,f,g,LSfunEvals] = ArmijoBacktrack(x,t,d,f,fr,g,gtd,c1,LS,tolX,debug,doPlot,1,funObj,varargin{:});
        end
        funEvals = funEvals + LSfunEvals;

    elseif LS < 6
        % Find Point satisfying Wolfe

        if computeHessian
            [t,f,g,LSfunEvals,H] = WolfeLineSearch(x,t,d,f,g,gtd,c1,c2,LS,25,tolX,debug,doPlot,LS_saveHessianComp,funObj,varargin{:});
        else
            [t,f,g,LSfunEvals] = WolfeLineSearch(x,t,d,f,g,gtd,c1,c2,LS,25,tolX,debug,doPlot,1,funObj,varargin{:});
        end
        funEvals = funEvals + LSfunEvals;
        x = x + t*d;

    else
        % Use Matlab optim toolbox line search
        [t,f_new,fPrime_new,g_new,LSexitFlag,LSiter]=...
            lineSearch({'fungrad',[],funObj},x,p,1,p,d,f,gtd,t,c1,c2,-inf,maxFunEvals-funEvals,...
            tolX,[],[],[],varargin{:});
        funEvals = funEvals + LSiter;
        if isempty(t)
            exitflag = -2;
            msg = 'Matlab LineSearch failed';
            break;
        end

        if method >= NEWTON
            [f_new,g_new,H] = funObj(x + t*d,varargin{:});
            funEvals = funEvals + 1;
        end
        x = x + t*d;
        f = f_new;
        g = g_new;
    end

    % Output iteration information
    if verboseI
        fprintf('%10d %10d %15.5e %15.5e %15.5e\n',i,funEvals*funEvalMultiplier,t,f,sum(abs(g)));
    end

    % Output Function
    if ~isempty(outputFcn)
        callOutput(outputFcn,x,'iter',i,funEvals,f,t,gtd,g,d,sum(abs(g)),varargin{:});
    end

    % Update Trace
    trace.fval(end+1,1) = f;
    trace.funcCount(end+1,1) = funEvals;

    % Check Optimality Condition
    if sum(abs(g)) <= tolFun
        exitflag=1;
        msg = 'Optimality Condition below TolFun';
        break;
    end

    % ******************* Check for lack of progress *******************

    if sum(abs(t*d)) <= tolX
        exitflag=2;
        msg = 'Step Size below TolX';
        break;
    end

    if abs(f-f_old) < tolX
        exitflag=2;
        msg = 'Function Value changing by less than TolX';
        break;
    end

    % ******** Check for going over iteration/evaluation limit *******************

    if funEvals*funEvalMultiplier > maxFunEvals
        exitflag = 0;
        msg = 'Exceeded Maximum Number of Function Evaluations';
        break;
    end

    if i == maxIter
        exitflag = 0;
        msg='Exceeded Maximum Number of Iterations';
        break;
    end

end

if verbose
    fprintf('%s\n',msg);
end
if nargout > 3
    output = struct('iterations',i,'funcCount',funEvals*funEvalMultiplier,...
        'algorithm',method,'firstorderopt',sum(abs(g)),'message',msg,'trace',trace);
end

% Output Function
if ~isempty(outputFcn)
    callOutput(outputFcn,x,'done',i,funEvals,f,t,gtd,g,d,sum(abs(g)),varargin{:});
end

end

author: huashiyiqike

repo: LSTM-MATLAB


text:
%MULTIPLY  Calculates the per-element scaled product of two arrays
%
%     dst = cv.multiply(src1, src2)
%     dst = cv.multiply(src1, src2, 'OptionName',optionValue, ...)
%
% ## Input
% * __src1__ first input array.
% * __src2__ second input array of the same size and type as `src1`.
%
% ## Output
% * __dst__ output array of the same size and type as `src1`.
%
% ## Options
% * __Scale__ optional scalar factor. default 1
% * __DType__ optional depth of the output array: `uint8`, `int16`, `double`,
%   etc. default -1
%
% The function cv.multiply calculates the per-element product of two arrays:
%
%     dst(I) = saturate(scale*src1(I) * src2(I))
%
% Note: Saturation is not applied when the output array has the depth `int32`.
% You may even get result of an incorrect sign in the case of overflow.
%
% See also: cv.add, cv.subtract, cv.divide, cv.addWeighted, cv.accumulate,
%  cv.accumulateProduct, cv.accumulateSquare, immultiply, times, mtimes
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% SaveCalPar   save calibration parameters in different formats
%
% [Cst,Rot] = savecalpar(P,config)
% P ... 3*CAM x 4 matrix containing result of selfcalibration, see EUCLIDIZE
% config ... configuration structure, see CONFIGDATA
%
%
% Cst ... CAMSx3   matrix containing the camera centers (in world coord.)
% Rot ... 3*CAMSx3 matrix containing camera rotation matrices

% $Author: svoboda $
% $Revision: 2.0 $
% $Id: savecalpar.m,v 2.0 2003/06/19 12:07:03 svoboda Exp $
% $State: Exp $

function [Cst,Rot] = savecalpar(P,config)

Octave = exist('OCTAVE_VERSION', 'builtin') ~= 0;

idxused = config.cal.cams2use;

CAMS = size(P,1)/3;

Cst = zeros(CAMS,3);
Pst = zeros(3*CAMS,3);
Rot = [];
for i=1:CAMS,
  % do not save P matrices in separate files
  if 1
	Pmat = P(i*3-2:i*3,:);
        if Octave
          save(sprintf(config.files.CalPmat,idxused(i)),'Pmat');
        else
          save(sprintf(config.files.CalPmat,idxused(i)),'Pmat','-ASCII');
        end
  end
  sc = norm(P(i*3,1:3));
  % first normalize the Projection matrices to get normalized pixel points
  P(i*3-2:i*3,:) = P(i*3-2:i*3,:)./sc;
  % decompose the matrix by using rq decomposition
  [K,R] = rq(P(i*3-2:i*3,1:3));
  tvec= inv(K)*P(i*3-2:i*3,4);			% translation vector
  C	  = -R'*tvec;						% camera center
  % Stephi calib params
  Pstephi		   = R'*inv(K);
  Pst(i*3-2:i*3,:) = Pstephi;
  Cst(i,:)		   = C';
  % Stephi requires to save the pars in more "wordy" form
  fid = fopen(sprintf(config.files.StCalPar,idxused(i)),'wt');
  if ~fid
	error('SaveCalPar: The camera cal file cannot be opened');
  else
	fprintf(fid,'C1 = %f \n', C(1));
	fprintf(fid,'C2 = %f \n', C(2));
	fprintf(fid,'C3 = %f \n', C(3));
	fprintf(fid,'\n');
	fprintf(fid,'P11 = %f \n', Pstephi(1,1));
	fprintf(fid,'P12 = %f \n', Pstephi(1,2));
	fprintf(fid,'P13 = %f \n', Pstephi(1,3));
	fprintf(fid,'P21 = %f \n', Pstephi(2,1));
	fprintf(fid,'P22 = %f \n', Pstephi(2,2));
	fprintf(fid,'P23 = %f \n', Pstephi(2,3));
	fprintf(fid,'P31 = %f \n', Pstephi(3,1));
	fprintf(fid,'P32 = %f \n', Pstephi(3,2));
	fprintf(fid,'P33 = %f \n', Pstephi(3,3));
	fclose(fid);
  end
  Rot	 = [Rot;R];
  % Prithwijit requires to save the pars in more "wordy" form
  if 0 % do not save in the Prithwijit format
	fid = fopen(sprintf(config.files.CalPar,idxused(i)),'wt');
	if ~fid
	  error('SaveCalPar: The camera cal file cannot be opened');
	else
	  fprintf(fid,'R11 = %f \n',R(1,1));
	  fprintf(fid,'R12 = %f \n',R(1,2));
	  fprintf(fid,'R13 = %f \n',R(1,3));
	  fprintf(fid,'R21 = %f \n',R(2,1));
	  fprintf(fid,'R22 = %f \n',R(2,2));
	  fprintf(fid,'R23 = %f \n',R(2,3));
	  fprintf(fid,'R31 = %f \n',R(3,1));
	  fprintf(fid,'R32 = %f \n',R(3,2));
	  fprintf(fid,'R33 = %f \n',R(3,3));
	  fprintf(fid,'\n');
	  fprintf(fid,'t11 = %f \n',tvec(1));
	  fprintf(fid,'t21 = %f \n',tvec(2));
	  fprintf(fid,'t31 = %f \n',tvec(3));
	  fprintf(fid,'\n');
	  fprintf(fid,'K11 = %f \n',K(1,1));
	  fprintf(fid,'K12 = %f \n',K(1,2));
	  fprintf(fid,'K13 = %f \n',K(1,3));
	  fprintf(fid,'K21 = %f \n',K(2,1));
	  fprintf(fid,'K22 = %f \n',K(2,2));
	  fprintf(fid,'K23 = %f \n',K(2,3));
	  fprintf(fid,'K31 = %f \n',K(3,1));
	  fprintf(fid,'K32 = %f \n',K(3,2));
	  fprintf(fid,'K33 = %f \n',K(3,3));
	  fprintf(fid,'\n');
	  fclose(fid);
	end
  end
end

% save Stehpi params
if Octave
  save(config.files.Pst,'Pst');
  save(config.files.Cst,'Cst');
else
  save(config.files.Pst,'Pst','-ASCII');
  save(config.files.Cst,'Cst','-ASCII');
end

author: strawlab

repo: MultiCamSelfCal


text:
function divdif_test18 ( )

%*****************************************************************************80
%
%% DIVDIF_TEST18 tests ROOTS_TO_DIF and DIF_TO_R8POLY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 October 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'DIVDIF_TEST18\n' );
  fprintf ( 1, '  ROOTS_TO_DIF computes the divided difference\n' );
  fprintf ( 1, '  polynomial with given roots;\n' );
  fprintf ( 1, '  DIF_TO_R8POLY converts it to a standard form\n' );
  fprintf ( 1, '  polynomial.\n' );
  fprintf ( 1, '\n' );

  nroots = 1;
  roots(1) = 3.0;
  r8vec_print ( nroots, roots, '  The roots:' );

  [ ntab, xtab, diftab ] = roots_to_dif ( nroots, roots );

  c = dif_to_r8poly ( ntab, xtab, diftab );

  r8poly_print ( ntab, c, '  The polynomial:' );

  nroots = 2;
  roots(1:2) = [ 3.0, 1.0 ];
  r8vec_print ( nroots, roots, '  The roots:' );

  [ ntab, xtab, diftab ] = roots_to_dif ( nroots, roots );

  c = dif_to_r8poly ( ntab, xtab, diftab );

  r8poly_print ( ntab, c, '  The polynomial:' );

  nroots = 3;
  roots(1:3) = [ 3.0, 1.0, 2.0 ];
  r8vec_print ( nroots, roots, '  The roots:' );

  [ ntab, xtab, diftab ] = roots_to_dif ( nroots, roots );

  c = dif_to_r8poly ( ntab, xtab, diftab );

  r8poly_print ( ntab, c, '  The polynomial:' );

  nroots = 4;
  roots(1:4) = [ 3.0, 1.0, 2.0, 4.0 ];
  r8vec_print ( nroots, roots, '  The roots:' );

  [ ntab, xtab, diftab ] = roots_to_dif ( nroots, roots );

  c = dif_to_r8poly ( ntab, xtab, diftab );

  r8poly_print ( ntab, c, '  The polynomial:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function class_info=gen_class_info_voc()

class_info=[];

class_info.class_names = { 'background', 'aeroplane', 'bicycle', 'bird', 'boat',  'bottle', 'bus',...
            'car', 'cat', 'chair', 'cow', 'diningtable','dog', 'horse',...
            'motorbike', 'person', 'pottedplant', 'sheep', 'sofa', 'train', 'tvmonitor', ...
            'void'}';

void_class_value=255;
class_info.class_label_values=uint8([0:20 void_class_value]);
class_info.background_label_value=uint8(0);
class_info.void_label_values=uint8(void_class_value);

% addpath ../libs/VOCdevkit_2012/VOCcode
class_info.mask_cmap = VOClabelcolormap(256);

class_info=process_class_info(class_info);

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function computeAllPrediction_batchLGG(pathExperiments,maxOrder,nBoot,imbalance,nBatch,matlabPATH,seed)
% -------------------------------------------------------------------------
% function computeAllPrediction_batchLGG(pathExperiments,maxOrder,nBoot,imbalance,nBatch,matlabPATH,seed)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes prediction performance estimation for a given
% feature set type, and for all model orders of all experiments with
% different degrees of freedom. See ref. [1] for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathExperiments: Full path to where all experiments need to be
%                     performed.
%                     --> Ex: /myProject/WORKSPACE/LOGISTIC_REGRESSION
% 2. maxOrder: Integer specifying the maximal model order to construct.
%              --> Ex: 10
% 3. nBoot: Number of bootstrap samples to use.
%           --> Ex: 100
% 4. imbalance: String specifying the type of imbalance-adjustement strategy
%               employed. Either 'IABR' for imbalance-adjusted bootstrap
%               resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%               logistic regression (formal reference to come).
%               --> Ex: 'IALR'
% 5. nBatch: Number of parallel batch.
%            --> Ex: 8
% 6. matlabPATH: Full path to the MATLAB excutable on the system.
%                --> 'matlab' if a symbolic link to the matlab executable
%                     was previously created.
% 7. seed: Random generator seed for reproducibility of experiments.
% -------------------------------------------------------------------------
% OUTPUTS: Prediction performance results are saved in a folder named 'RESULTS' in the
% corresponding folder of 'pathExperiments'.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;

% INITIALIZATON
time = 60; % Number of seconds to wait before checking if parallel computations are done
cd(pathExperiments), load('training')
pathModels = fullfile(pwd,'MODELS'); mkdir('RESULTS'), cd('RESULTS'), pathResults = pwd;
mkdir('batchLog_Results'), cd('batchLog_Results'), pathBatch = pwd;
nameOutcomes = fieldnames(training); nOutcomes = numel(nameOutcomes);
for o = 1:nOutcomes
    outcomes.(nameOutcomes{o}) = training.(nameOutcomes{o}).outcome;
end
setNames = fieldnames(training.(nameOutcomes{1}).text);
[param] = batchExperiments(setNames,outcomes,nBatch); nBatch = length(param);

% PRODUCE BATCH COMPUTATIONS
save('workspace','pathModels','pathResults','training','param','maxOrder','nBoot','imbalance','seed'), pause(5);
for i = 1:nBatch
    nameScript = ['batch',num2str(i),'_script.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'load(''workspace'')\n');
    for j = 1:numel(param{i})
        fprintf(fid,['computeAllPrediction_LGG(pathModels,pathResults,training,param{',num2str(i),'}{',num2str(j),'},maxOrder,nBoot,imbalance,seed)\n']);
    end
    fprintf(fid,['system(''touch batch',num2str(i),'_end'');\n']);
    fprintf(fid,'clear all');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end

% WAITING LOOP
waitBatch(pathBatch,time,nBatch)
delete('workspace.mat')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [im_th,TOAref,trgt,ijdim_ref,ul,bbox,zen,azi,zc,B1Satu,B2Satu,B3Satu,resolu]=nd2toarbt(path_top,filename)
% convert DNs to TOA ref and BT
% Revisions:
% Add Landsat 9 (Shi 2/17/2022)
% Cannot use the GRIDobj to read our band, because it may result in Nan data.
% (Shi 2/26/2018)
% Added a target image based on blue band, which will be helpfull to
% resample and reproject the auxiliary data.  (Shi 2/16/2018)
% Use REF vs. DN instead of RAD vs. DN (Zhe 06/20/2013)
% Combined the Earth-Sun distance table into the function (Zhe 04/09/2013)
% Process Landsat 8 DN values (Zhe 04/04/2013)
% Proces the new metadata for Landsat TM/ETM+ images (Zhe 09/28/2012)
% Fixed bugs caused by earth-sun distance table (Zhe 01/15/2012)
%
% [im_th,TOAref,ijdim_ref,ul,zen,azi,zc,B1Satu,B2Satu,B3Satu,resolu]=nd2toarbt(filename)
% Where:
% Inputs:
% filename='L*MTL.txt';
% Outputs:
% 1) im_th = Brightness Temperature (BT)
% 2) TOAref = Top Of Atmoshpere (TOA) reflectance
% 3) trgt = Target Grid Object
% 4) ijdim = [nrows,ncols]; % dimension of optical bands
% 5) ul = [upperleft_mapx upperleft_mapy];
% 6) zen = solar zenith angle (degrees);
% 7) azi = solar azimuth angle (degrees);
% 8) zc = Zone Number
% 9,10,11) Saturation (true) in the Visible bands
% 12) resolution of Fmask results

% % % fprintf('Read in header information & TIF images\n');
[Lmax,Lmin,Qcalmax,Qcalmin,Refmax,Refmin,ijdim_ref,ijdim_thm,reso_ref,...
    reso_thm,ul,bbox,zen,azi,zc,Lnum,doy]=lndhdrread(path_top,filename);

% earth-sun distance see G. Chander et al. RSE 113 (2009) 893-903
Tab_ES_Dist = [
    1 0.98331
    2 0.98330
    3 0.98330
    4 0.98330
    5 0.98330
    6 0.98332
    7 0.98333
    8 0.98335
    9 0.98338
    10 0.98341
    11 0.98345
    12 0.98349
    13 0.98354
    14 0.98359
    15 0.98365
    16 0.98371
    17 0.98378
    18 0.98385
    19 0.98393
    20 0.98401
    21 0.98410
    22 0.98419
    23 0.98428
    24 0.98439
    25 0.98449
    26 0.98460
    27 0.98472
    28 0.98484
    29 0.98496
    30 0.98509
    31 0.98523
    32 0.98536
    33 0.98551
    34 0.98565
    35 0.98580
    36 0.98596
    37 0.98612
    38 0.98628
    39 0.98645
    40 0.98662
    41 0.98680
    42 0.98698
    43 0.98717
    44 0.98735
    45 0.98755
    46 0.98774
    47 0.98794
    48 0.98814
    49 0.98835
    50 0.98856
    51 0.98877
    52 0.98899
    53 0.98921
    54 0.98944
    55 0.98966
    56 0.98989
    57 0.99012
    58 0.99036
    59 0.99060
    60 0.99084
    61 0.99108
    62 0.99133
    63 0.99158
    64 0.99183
    65 0.99208
    66 0.99234
    67 0.99260
    68 0.99286
    69 0.99312
    70 0.99339
    71 0.99365
    72 0.99392
    73 0.99419
    74 0.99446
    75 0.99474
    76 0.99501
    77 0.99529
    78 0.99556
    79 0.99584
    80 0.99612
    81 0.99640
    82 0.99669
    83 0.99697
    84 0.99725
    85 0.99754
    86 0.99782
    87 0.99811
    88 0.99840
    89 0.99868
    90 0.99897
    91 0.99926
    92 0.99954
    93 0.99983
    94 1.00012
    95 1.00041
    96 1.00069
    97 1.00098
    98 1.00127
    99 1.00155
    100 1.00184
    101 1.00212
    102 1.00240
    103 1.00269
    104 1.00297
    105 1.00325
    106 1.00353
    107 1.00381
    108 1.00409
    109 1.00437
    110 1.00464
    111 1.00492
    112 1.00519
    113 1.00546
    114 1.00573
    115 1.00600
    116 1.00626
    117 1.00653
    118 1.00679
    119 1.00705
    120 1.00731
    121 1.00756
    122 1.00781
    123 1.00806
    124 1.00831
    125 1.00856
    126 1.00880
    127 1.00904
    128 1.00928
    129 1.00952
    130 1.00975
    131 1.00998
    132 1.01020
    133 1.01043
    134 1.01065
    135 1.01087
    136 1.01108
    137 1.01129
    138 1.01150
    139 1.01170
    140 1.01191
    141 1.01210
    142 1.01230
    143 1.01249
    144 1.01267
    145 1.01286
    146 1.01304
    147 1.01321
    148 1.01338
    149 1.01355
    150 1.01371
    151 1.01387
    152 1.01403
    153 1.01418
    154 1.01433
    155 1.01447
    156 1.01461
    157 1.01475
    158 1.01488
    159 1.01500
    160 1.01513
    161 1.01524
    162 1.01536
    163 1.01547
    164 1.01557
    165 1.01567
    166 1.01577
    167 1.01586
    168 1.01595
    169 1.01603
    170 1.01610
    171 1.01618
    172 1.01625
    173 1.01631
    174 1.01637
    175 1.01642
    176 1.01647
    177 1.01652
    178 1.01656
    179 1.01659
    180 1.01662
    181 1.01665
    182 1.01667
    183 1.01668
    184 1.01670
    185 1.01670
    186 1.01670
    187 1.01670
    188 1.01669
    189 1.01668
    190 1.01666
    191 1.01664
    192 1.01661
    193 1.01658
    194 1.01655
    195 1.01650
    196 1.01646
    197 1.01641
    198 1.01635
    199 1.01629
    200 1.01623
    201 1.01616
    202 1.01609
    203 1.01601
    204 1.01592
    205 1.01584
    206 1.01575
    207 1.01565
    208 1.01555
    209 1.01544
    210 1.01533
    211 1.01522
    212 1.01510
    213 1.01497
    214 1.01485
    215 1.01471
    216 1.01458
    217 1.01444
    218 1.01429
    219 1.01414
    220 1.01399
    221 1.01383
    222 1.01367
    223 1.01351
    224 1.01334
    225 1.01317
    226 1.01299
    227 1.01281
    228 1.01263
    229 1.01244
    230 1.01225
    231 1.01205
    232 1.01186
    233 1.01165
    234 1.01145
    235 1.01124
    236 1.01103
    237 1.01081
    238 1.01060
    239 1.01037
    240 1.01015
    241 1.00992
    242 1.00969
    243 1.00946
    244 1.00922
    245 1.00898
    246 1.00874
    247 1.00850
    248 1.00825
    249 1.00800
    250 1.00775
    251 1.00750
    252 1.00724
    253 1.00698
    254 1.00672
    255 1.00646
    256 1.00620
    257 1.00593
    258 1.00566
    259 1.00539
    260 1.00512
    261 1.00485
    262 1.00457
    263 1.00430
    264 1.00402
    265 1.00374
    266 1.00346
    267 1.00318
    268 1.00290
    269 1.00262
    270 1.00234
    271 1.00205
    272 1.00177
    273 1.00148
    274 1.00119
    275 1.00091
    276 1.00062
    277 1.00033
    278 1.00005
    279 0.99976
    280 0.99947
    281 0.99918
    282 0.99890
    283 0.99861
    284 0.99832
    285 0.99804
    286 0.99775
    287 0.99747
    288 0.99718
    289 0.99690
    290 0.99662
    291 0.99634
    292 0.99605
    293 0.99577
    294 0.99550
    295 0.99522
    296 0.99494
    297 0.99467
    298 0.99440
    299 0.99412
    300 0.99385
    301 0.99359
    302 0.99332
    303 0.99306
    304 0.99279
    305 0.99253
    306 0.99228
    307 0.99202
    308 0.99177
    309 0.99152
    310 0.99127
    311 0.99102
    312 0.99078
    313 0.99054
    314 0.99030
    315 0.99007
    316 0.98983
    317 0.98961
    318 0.98938
    319 0.98916
    320 0.98894
    321 0.98872
    322 0.98851
    323 0.98830
    324 0.98809
    325 0.98789
    326 0.98769
    327 0.98750
    328 0.98731
    329 0.98712
    330 0.98694
    331 0.98676
    332 0.98658
    333 0.98641
    334 0.98624
    335 0.98608
    336 0.98592
    337 0.98577
    338 0.98562
    339 0.98547
    340 0.98533
    341 0.98519
    342 0.98506
    343 0.98493
    344 0.98481
    345 0.98469
    346 0.98457
    347 0.98446
    348 0.98436
    349 0.98426
    350 0.98416
    351 0.98407
    352 0.98399
    353 0.98391
    354 0.98383
    355 0.98376
    356 0.98370
    357 0.98363
    358 0.98358
    359 0.98353
    360 0.98348
    361 0.98344
    362 0.98340
    363 0.98337
    364 0.98335
    365 0.98333
    366 0.98331];

if Lnum >= 4 && Lnum <= 7
    % LPGS Upper lef corner alignment (see Landsat handbook for detail)
    ul(1)=ul(1)-15;
    ul(2)=ul(2)+15;
    resolu=[reso_ref,reso_ref];
    % Read in all bands
    % Band1
    n_B1=dir(fullfile(path_top,'L*B1*'));
    if isempty(n_B1)
        n_B1=dir(fullfile(path_top,'L*b1*'));
    end
    im_B1=single(imread(fullfile(path_top,n_B1.name)));
    % Band2
    % Also used to be target or referenced image when resampling and reprojecting to the image extent and resolution.
    n_B2=dir(fullfile(path_top,'L*B2*'));
    if isempty(n_B2)
        n_B2=dir(fullfile(path_top,'L*b2*'));
    end
    im_B2=single(imread(fullfile(path_top,n_B2.name)));
    trgt = GRIDobj(fullfile(path_top,n_B2.name));
%     im_B2=single(trgt.Z);
    trgt.Z=[];% remove non-used values to save memory.

    % Band3
    n_B3=dir(fullfile(path_top,'L*B3*'));
    if isempty(n_B3)
        n_B3=dir(fullfile(path_top,'L*b3*'));
    end
    im_B3=single(imread(fullfile(path_top,n_B3.name)));
    % Band4
    n_B4=dir(fullfile(path_top,'L*B4*'));
    if isempty(n_B4)
        n_B4=dir(fullfile(path_top,'L*b4*'));
    end
    im_B4=single(imread(fullfile(path_top,n_B4.name)));
    % Band5
    n_B5=dir(fullfile(path_top,'L*B5*'));
    if isempty(n_B5)
        n_B5=dir(fullfile(path_top,'L*b5*'));
    end
    im_B5=single(imread(fullfile(path_top,n_B5.name)));
    % Band6
    if Lnum==7
        n_B6=dir(fullfile(path_top,'L*B6*1*'));
        if isempty(n_B6)
            n_B6=dir(fullfile(path_top,'L*b6*1*'));
        end
    else
        n_B6=dir(fullfile(path_top,'L*B6*'));
        if isempty(n_B6)
            n_B6=dir(fullfile(path_top,'L*b6*'));
        end
    end
    im_th=single(imread(fullfile(path_top,n_B6.name)));
    % check to see whether need to resample thermal band
    if reso_ref~=reso_thm
        % resmaple thermal band
        im_th=pixel2pixv([ul(2),ul(1)],[ul(2),ul(1)],...
            resolu,[reso_thm,reso_thm],...
            im_th,[ijdim_ref(2),ijdim_ref(1)],[ijdim_thm(2),ijdim_thm(1)]);
    end
    % Band7
    n_B7=dir(fullfile(path_top,'L*B7*'));
    if isempty(n_B7)
        n_B7=dir(fullfile(path_top,'L*b7*'));
    end
    im_B7=single(imread(fullfile(path_top,n_B7.name)));
    % only processing pixesl where all bands have values (id_mssing)
    id_missing=im_B1==0|im_B2==0|im_B3==0|im_B4==0|im_B5==0|im_th==0|im_B7==0;
    % find pixels that are saturated in the visible bands
    B1Satu=im_B1==255;
    B2Satu=im_B2==255;
    B3Satu=im_B3==255;

    % ND to radiance first
% % %     fprintf('From DNs to TOA ref & BT\n');
    im_B1=((Lmax(1)-Lmin(1))/(Qcalmax(1)-Qcalmin(1)))*(im_B1-Qcalmin(1))+Lmin(1);
    im_B2=((Lmax(2)-Lmin(2))/(Qcalmax(2)-Qcalmin(2)))*(im_B2-Qcalmin(2))+Lmin(2);
    im_B3=((Lmax(3)-Lmin(3))/(Qcalmax(3)-Qcalmin(3)))*(im_B3-Qcalmin(3))+Lmin(3);
    im_B4=((Lmax(4)-Lmin(4))/(Qcalmax(4)-Qcalmin(4)))*(im_B4-Qcalmin(4))+Lmin(4);
    im_B5=((Lmax(5)-Lmin(5))/(Qcalmax(5)-Qcalmin(5)))*(im_B5-Qcalmin(5))+Lmin(5);
    im_th=((Lmax(6)-Lmin(6))/(Qcalmax(6)-Qcalmin(6)))*(im_th-Qcalmin(6))+Lmin(6);
    im_B7=((Lmax(7)-Lmin(7))/(Qcalmax(7)-Qcalmin(7)))*(im_B7-Qcalmin(7))+Lmin(7);

    % Landsat 4, 5, and 7 solar spectral irradiance
    % see G. Chander et al. RSE 113 (2009) 893-903
    esun_L7=[1997.000, 1812.000, 1533.000, 1039.000, 230.800, -1.0, 84.90];
    esun_L5=[1983.0, 1796.0, 1536.0, 1031.0, 220.0, -1.0, 83.44];
    esun_L4=[1983.0, 1795.0, 1539.0, 1028.0, 219.8, -1.0, 83.49];

    if Lnum==7
        ESUN=esun_L7;
    elseif Lnum==5
        ESUN=esun_L5;
    elseif Lnum==4
        ESUN=esun_L4;
    end

    % earth-sun distance see G. Chander et al. RSE 113 (2009) 893-903
    dsun_doy = Tab_ES_Dist(doy,2);

    % compute TOA reflectances
    % converted from degrees to radiance
    s_zen=deg2rad(zen);
    im_B1=10^4*pi*im_B1*dsun_doy^2/(ESUN(1)*cos(s_zen));
    im_B2=10^4*pi*im_B2*dsun_doy^2/(ESUN(2)*cos(s_zen));
    im_B3=10^4*pi*im_B3*dsun_doy^2/(ESUN(3)*cos(s_zen));
    im_B4=10^4*pi*im_B4*dsun_doy^2/(ESUN(4)*cos(s_zen));
    im_B5=10^4*pi*im_B5*dsun_doy^2/(ESUN(5)*cos(s_zen));
    im_B7=10^4*pi*im_B7*dsun_doy^2/(ESUN(7)*cos(s_zen));

    % convert Band6 from radiance to BT
    % fprintf('From Band 6 Radiance to Brightness Temperature\n');
    % see G. Chander et al. RSE 113 (2009) 893-903
    K1_L4=  671.62;
    K2_L4= 1284.30;
    K1_L5=  607.76;
    K2_L5= 1260.56;
    K1_L7 =  666.09;
    K2_L7 = 1282.71;

    if Lnum==7
        K1=K1_L7;
        K2=K2_L7;
    elseif Lnum==5
        K1=K1_L5;
        K2=K2_L5;
    elseif Lnum==4
        K1=K1_L4;
        K2=K2_L4;
    end

    im_th=K2./log((K1./im_th)+1);
    % convert from Kelvin to Celcius with 0.01 scale_facor
    im_th=100*(im_th-273.15);
    % get data ready for Fmask
    im_B1(id_missing)=-9999;
    im_B2(id_missing)=-9999;
    im_B3(id_missing)=-9999;
    im_B4(id_missing)=-9999;
    im_B5(id_missing)=-9999;
    im_th(id_missing)=-9999;
    im_B7(id_missing)=-9999;
    clear id_missing;
    TOAref=zeros(ijdim_ref(1),ijdim_ref(2),6,'single');% Band1,2,3,4,5,&7
    TOAref(:,:,1)=im_B1;clear im_B1;
    TOAref(:,:,2)=im_B2;clear im_B2;
    TOAref(:,:,3)=im_B3;clear im_B3;
    TOAref(:,:,4)=im_B4;clear im_B4;
    TOAref(:,:,5)=im_B5;clear im_B5;
    TOAref(:,:,6)=im_B7;clear im_B7;

elseif Lnum == 8 || Lnum == 9
    % LPGS Upper lef corner alignment (see Landsat handbook for detail)
    ul(1)=ul(1)-15;
    ul(2)=ul(2)+15;
    resolu=[reso_ref,reso_ref];
    % Read in all bands
    % Band2
    n_B2=dir(fullfile(path_top,'L*B2*'));
    if isempty(n_B2)
        n_B2=dir(fullfile(path_top,'L*b2*'));
    end
    im_B2=single(imread(fullfile(path_top,n_B2(1).name))); clear n_B2;
    % Band3
    % Also used to be target or referenced image when resampling and reprojecting to the image extent and resolution.
    n_B3=dir(fullfile(path_top,'L*B3*'));
    if isempty(n_B3)
        n_B3=dir(fullfile(path_top,'L*b3*'));
    end
    im_B3=single(imread(fullfile(path_top,n_B3(1).name)));
    trgt = GRIDobj(fullfile(path_top,n_B3.name)); clear n_B3;
%     im_B3=single(trgt.Z);
    trgt.Z=[];% remove non-used values to save memory.

    % Band4
    n_B4=dir(fullfile(path_top,'L*B4*'));
    if isempty(n_B4)
        n_B4=dir(fullfile(path_top,'L*b4*'));
    end
    im_B4=single(imread(fullfile(path_top,n_B4(1).name))); clear n_B4;
    % Band5
    n_B5=dir(fullfile(path_top,'L*B5*'));
    if isempty(n_B5)
        n_B5=dir(fullfile(path_top,'L*b5*'));
    end
    im_B5=single(imread(fullfile(path_top,n_B5(1).name))); clear n_B5;
    % Band6
    n_B6=dir(fullfile(path_top,'L*B6*'));
    if isempty(n_B6)
        n_B6=dir(fullfile(path_top,'L*b6*'));
    end
    im_B6=single(imread(fullfile(path_top,n_B6(1).name))); clear n_B6;
    % Band7
    n_B7=dir(fullfile(path_top,'L*B7*'));
    if isempty(n_B7)
        n_B7=dir(fullfile(path_top,'L*b7*'));
    end
    im_B7=single(imread(fullfile(path_top,n_B7(1).name))); clear n_B7;
    % Band9
    n_B9=dir(fullfile(path_top,'L*B9.*'));
    if isempty(n_B9)
        n_B9=dir(fullfile(path_top,'L*b9.*'));
    end
    im_B9=single(imread(fullfile(path_top,n_B9(1).name))); clear n_B9;
    % Band10
    n_B10=dir(fullfile(path_top,'L*B10*'));
    if isempty(n_B10)
        n_B10=dir(fullfile(path_top,'L*b10*'));
    end
    im_B10=single(imread(fullfile(path_top,n_B10(1).name))); clear n_B10 path_top;

    % check to see whether need to resample thermal band
    if reso_ref~=reso_thm
        % resmaple thermal band
        im_B10=pixel2pixv([ul(2),ul(1)],[ul(2),ul(1)],...
            resolu,[reso_thm,reso_thm],...
            im_B10,[ijdim_ref(2),ijdim_ref(1)],[ijdim_thm(2),ijdim_thm(1)]);
    end

    % only processing pixesl where all bands have values (id_mssing)
    id_missing=im_B2==0|im_B3==0|im_B4==0|im_B5==0|im_B6==0|im_B7==0|im_B9==0|im_B10==0;
    % find pixels that are saturated in the visible bands
    B1Satu=im_B2==65535;
    B2Satu=im_B3==65535;
    B3Satu=im_B4==65535;

    % ND to TOA reflectance with 0.0001 scale_facor
% % %     fprintf('From DNs to TOA ref & BT\n');
    im_B2=((Refmax(1)-Refmin(1))/(Qcalmax(1)-Qcalmin(1)))*(im_B2-Qcalmin(1))+Refmin(1);
    im_B3=((Refmax(2)-Refmin(2))/(Qcalmax(2)-Qcalmin(2)))*(im_B3-Qcalmin(2))+Refmin(2);
    im_B4=((Refmax(3)-Refmin(3))/(Qcalmax(3)-Qcalmin(3)))*(im_B4-Qcalmin(3))+Refmin(3);
    im_B5=((Refmax(4)-Refmin(4))/(Qcalmax(4)-Qcalmin(4)))*(im_B5-Qcalmin(4))+Refmin(4);
    im_B6=((Refmax(5)-Refmin(5))/(Qcalmax(5)-Qcalmin(5)))*(im_B6-Qcalmin(5))+Refmin(5);
    im_B7=((Refmax(6)-Refmin(6))/(Qcalmax(6)-Qcalmin(6)))*(im_B7-Qcalmin(6))+Refmin(6);
    im_B9=((Refmax(7)-Refmin(7))/(Qcalmax(7)-Qcalmin(7)))*(im_B9-Qcalmin(7))+Refmin(7);
    im_B10=((Lmax(8)-Lmin(8))/(Qcalmax(8)-Qcalmin(8)))*(im_B10-Qcalmin(8))+Lmin(8);

    % compute TOA reflectances
    % with a correction for the sun angle
    s_zen=deg2rad(zen);
    im_B2=10^4*im_B2/cos(s_zen);
    im_B3=10^4*im_B3/cos(s_zen);
    im_B4=10^4*im_B4/cos(s_zen);
    im_B5=10^4*im_B5/cos(s_zen);
    im_B6=10^4*im_B6/cos(s_zen);
    im_B7=10^4*im_B7/cos(s_zen);
    im_B9=10^4*im_B9/cos(s_zen);

    % convert Band6 from radiance to BT
    % fprintf('From Band 6 Radiance to Brightness Temperature\n');
    K1_B10 =  774.89;
    K2_B10 = 1321.08;

    im_B10=K2_B10./log((K1_B10./im_B10)+1);

    % convert from Kelvin to Celcius with 0.01 scale_facor
    im_B10=100*(im_B10-273.15);

    % get data ready for Fmask
    im_B2(id_missing)=-9999;
    im_B3(id_missing)=-9999;
    im_B4(id_missing)=-9999;
    im_B5(id_missing)=-9999;
    im_B6(id_missing)=-9999;
    im_B7(id_missing)=-9999;
    im_B9(id_missing)=-9999;
    im_B10(id_missing) = -9999;
    clear id_missing; % empty memory

    TOAref=zeros(ijdim_ref(1),ijdim_ref(2),7,'single');% Band 2,3,4,5,6,7,& 9
    TOAref(:,:,1)=im_B2; clear im_B2;
    TOAref(:,:,2)=im_B3; clear im_B3;
    TOAref(:,:,3)=im_B4; clear im_B4;
    TOAref(:,:,4)=im_B5; clear im_B5;
    TOAref(:,:,5)=im_B6; clear im_B6;
    TOAref(:,:,6)=im_B7; clear im_B7;
    TOAref(:,:,7)=im_B9; clear im_B9;

    im_th=zeros(ijdim_ref(1),ijdim_ref(2),'single');% Band 10
    im_th(:,:) = im_B10; clear im_B10;
else
    fprintf('This sensor is not Landsat 4, 5, 7, or 8!\n');
    return;
end

end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
classdef VaEA < ALGORITHM
% <multi/many> <real/integer/label/binary/permutation>
% Vector angle based evolutionary algorithm

%------------------------------- Reference --------------------------------
% Y. Xiang, Y. Zhou, M. Li, and Z. Chen, A vector angle-based evolutionary
% algorithm for unconstrained many-objective optimization, IEEE
% Transactions on Evolutionary Computation, 2017, 21(1): 131-152.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = randi(Problem.N,1,Problem.N);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                Population = EnvironmentalSelection([Population,Offspring],Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function modifieddq = satdq(dqlimit,q)
modifieddq = q;
count = length(q);
for i=1:count
    if modifieddq(i)<-dqlimit(i)
        modifieddq(i) = -dqlimit(i);
    end
    if modifieddq(i)> dqlimit(i)
        modifieddq(i) = dqlimit(i);
    end
end
end

author: xuhuairuogu

repo: V-REP-Simulation-Projects

path: github-repos/MATLAB/xuhuairuogu-V-REP-Simulation-Projects/V-REP-Simulation-Projects-841b944af4ea3a8fb250578d36434515f577f411/admittance control(adapted from an example in the book--A Systematic Approach to Learning Robot Programming with ROS)/iiwa14_kinematics/satdq.m


text:
function density = p12_density ( m, n, point )

%*****************************************************************************80
%
%% P12_DENSITY returns the density for problem 12.
%
%  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 N, the number of points.
%
%    Input, real POINT(M,N), the coordinates of the points.
%
%    Output, real DENSITY(N), the mesh density at each point.
%
  density = ones ( 1, n );

  return
end


text:
function [vert, faces] = tess_disc(h0)
% TESS_CIRCLE: Create a meshed disc.
%
% USAGE:  [vert, faces] = tess_disc(h0);
%
% INPUTS:
%    - h0 : Initial edge length

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2016

% Parse inputs
if (nargin < 1) || isempty(h0)
    h0 = 0.07;
end

% Distance and edge length function
fd = @(p) sqrt(sum(p.^2,2))-1;
fh = @(p) ones(size(p,1),1);
bbox = [-1,-1;1,1];
pfix = [];

% Use distmesh toolbox
[vert, faces] = distmesh2d(fd, fh, h0, bbox, pfix);

author: brainstorm-tools


text:
function d = dunion ( d1, d2 )

%*****************************************************************************80
%
%% DUNION returns the signed distance to a union of two regions.
%
%  Copyright:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  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, real D1, D2, the signed distance of one or more points
%    to two regions.
%
%    Output, real D, the signed distance of one or more points to
%    the region formed by the union of the two regions.
%
  d = min ( d1, d2 );

  return
end


text:
% Plot the curves extracted via curve_ext or curve_ext_multi
%
%  hc = plot_ext_curves(t, x, Tx, fs, Cs, Es, opt, clwin)
%
% Input:
%  t, x, opt: same as input to synsq_cwt_fw/iw
%  Tx, fs, clwin: same as input to curve_ext_multi
%  Cs, Es: same as output of curve_ext_multi
% Output:
%  hc: Object handle for the resulting figure
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------
function hc = plot_ext_curves(t, x, Tx, fs, Cs, Es, opt, clwin)

if nargin<8, clwin = 4; end
if nargin<7, opt = struct(); end

if ~isfield(opt, 'markers')
  opt.markers = {'+','o','*','x','s','d','^','v','>','<','p','h' };
end

Nc = length(Es);
[na, N] = size(Tx);

assert(N == length(t));
assert(N == length(x));
Nc = length(Es);

% Reconstruct the curve signals
xrs = curve_ext_recon(Tx, fs, Cs, opt, clwin);

hca = gcf;

% Plot everything
if isfield(opt, 'hc'),
    hc(1) = opt.hc(1);
    subplot(hc(1));
else
    hc(1) = subplot(2, 1, 1);
end

% Plot Tx
tplot(Tx, t, fs, opt);
% Set T-axis properly
axis tight;
xlabel('t');
ylabel('f');
title('Extracted Contour(s)');

% Plot the contours
hold on;
cpl=plot_markers(t, log2(fs(Cs)), opt.markers, '--k', 'LineWidth', 2);
legend(cpl, cellfun(@(j) {sprintf('k=%g', j)}, num2cell(1:Nc)));

if isfield(opt, 'hc')
    hc(2) = opt.hc(2);
    subplot(hc(2));
else
    hc(2) = subplot(2, 1, 2);
end

% % Plot both original signal and extracted contour
lcs = plot(t, x, 'k'); hold on;
lcs = [lcs plot_markers(t, xrs, opt.markers, '--k', 'LineWidth', 1)];

axis tight;
xlabel('t');
ylabel('x_c(t)');
grid on;
title('Reconstruction of Extracted Contour');
xcleg = arrayfun(@(j) {sprintf('x_%g(t)', j)}, 1:Nc);
legend(lcs, 'x(t)', xcleg{:});

linkaxes([hc(1) hc(2)], 'x');
zec = zoom(hca);
set(zec, 'RightClickAction', 'InverseZoom');
setAxesZoomMotion(zec, hc(1), 'vertical');
setAxesZoomMotion(zec, hc(2), 'horizontal');
set(zec, 'Enable', 'on');

end

author: HeLiangHIT


text:
function sparse_grid_mixed_size_tabulate ( rule_1d, alpha_1d, beta_1d, ...
   dim_min, dim_max, level_max_min, level_max_max )

%****************************************************************************80
%
%% SPARSE_GRID_MIXED_SIZE_TABULATE tests SPARSE_GRID_MIXED_SIZE.
%
%  Discussion:
%
%    We do NOT consider mixed rules.  Instead, we are looking at sparse grid
%    rules for which all dimensions use the same 1D rule family.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer RULE_1D, the 1D rule.
%
%    Input, real ALPHA_1D, BETA_1D, the optional parameters.
%
%    Input, integer DIM_MIN, the minimum spatial dimension.
%
%    Input, integer DIM_MAX, the maximum spatial dimension.
%
%    Input, integer LEVEL_MAX_MIN, the minimum value of LEVEL_MAX.
%
%    Input, integer LEVEL_MAX_MAX, the maximum value of LEVEL_MAX.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_MIXED_SIZE_TABULATE\n' );
  fprintf ( 1, '  SPARSE_GRID_MIXED_SIZE returns the number of distinct\n' );
  fprintf ( 1, '  points in a sparse grid.\n' );

  if ( rule_1d == 0 )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Here we report the total number of polynomials of\n' )
    fprintf ( 1, '  degree DEGREE or less in the given DIM dimensional space.\n' );
    fprintf ( 1, '  (This is essentially Pascal''s triangle.)\n' );
    fprintf ( 1, '\n' );

    fprintf ( 1, '   DIM: ' );
    for dim_num = dim_min : dim_max
      fprintf ( 1, '  %8d', dim_num );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '      DEGREE\n' );
    fprintf ( 1, '\n' );

    degree_max = 2 * level_max_max + 1;

    for degree = 0 : degree_max

      fprintf ( 1, '    %4d', degree );

      for dim_num = dim_min : dim_max

        point_num = i4_choose ( dim_num + degree, dim_num );

        fprintf ( 1, '  %8d', point_num );

      end

      fprintf ( 1, '\n' );

    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, '  We use the same rule in all dimensions, and count the points\n' );
    fprintf ( 1, '  for a range of dimensions and levels.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  1D rule index is %d\n', rule_1d );
    fprintf ( 1, '  ALPHA parameter is %f\n', alpha_1d );
    fprintf ( 1, '  BETA parameter is  %f\n', beta_1d );
    fprintf ( 1, '\n' );

    tol = sqrt ( eps );

    fprintf ( 1, '   DIM: ' );
    for dim_num = dim_min : dim_max
      fprintf ( 1, '  %8d', dim_num );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '   LEVEL_MAX\n' );
    fprintf ( 1, '\n' );

    for level_max = level_max_min : level_max_max

      fprintf ( 1, '    %4d', level_max );

      for dim_num = dim_min : dim_max

        rule(1:dim_num) = rule_1d;
        alpha(1:dim_num) = alpha_1d;
        beta(1:dim_num) = beta_1d;

        point_num = sparse_grid_mixed_size ( dim_num, level_max, rule, alpha, ...
          beta, tol );

        fprintf ( 1, '  %8d', point_num );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
function [ y, m, d, ierror ] = ymd_check_islamic ( y, m, d )

%*****************************************************************************80
%
%% YMD_CHECK_ISLAMIC checks an Islamic YMD date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, the YMD date, which may
%    be corrected if necessary and possible.
%
%    Output, integer IERROR, is 0 if the date is legal.
%

%
%  Check the year.
%
  if ( y <= 0 )
    ierror = 1;
    return
  end
%
%  Check the month.
%
  [ y, m, ierror ] = ym_check_islamic ( y, m );

  if ( ierror ~= 0 )
    return
  end
%
%  Check the day.
%
  [ y, m, d ] = day_borrow_islamic ( y, m, d );

  [ y, m, d ] = day_carry_islamic ( y, m, d );

  return
end


text:
function this = MrDataNd_constructor(this, testVariantsDataNd)
% Unit test for MrDataNd constructor with different input variants
%
%   Y = MrUnitTest()
%   run(Y, 'MrDataNd_constructor');
%
% This is a method of class MrUnitTest.
%
% IN
%
% OUT
%
% EXAMPLE
%   MrDataNd_constructor
%
%   See also MrUnitTest

% Author:   Saskia Bollmann
% Created:  2018-02-08
% Copyright (C) 2018 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.
%

% get the dimInfo reference
dimInfo = this.make_dimInfo_reference(0);
% change some parameters to check defaults are overwritten
dimInfo.dimLabels = {'dL1', 'dL2', 'dL3', 'dL4', 'dL5'};
dimInfo.units = {'u1', 'u2', 'u3', 'u4', 'u5'};
% remove samples to make data matrix smaller
dimInfo.dL4.nSamples = 10;
rng(1);
dataMatrix = rand(dimInfo.nSamples);

switch testVariantsDataNd
    case 'matrix'
        % load matrix from workspace
        actSolution = MrDataNd(dataMatrix);
        expSolution = actSolution;

    case 'matrixWithDimInfo'
        % load matrix with dimInfo object
        actSolution = MrDataNd(dataMatrix, 'dimInfo', dimInfo);
        expSolution = actSolution;
    case 'matrixWithPropVal'
        % load matrix with prop/val pairs describing dimInfo
        actSolution = MrDataNd(dataMatrix, 'dimLabels', dimInfo.dimLabels, ...
            'units', dimInfo.units, 'samplingPoints', dimInfo.samplingPoints);
        expSolution = actSolution;
end
% verify whether expected and actual solution are identical
% Note: convert dimInfo to struct, since the PublicPropertyComparator (to allow
% nans to be treated as equal) does not compare properties of objects that
% overload subsref

warning('off', 'MATLAB:structOnObject');
actSolution.dimInfo = struct(actSolution.dimInfo);
warning('on', 'MATLAB:structOnObject');
% verify equality of expected and actual solution
% import matlab.unittests to apply tolerances for objects
import matlab.unittest.TestCase
import matlab.unittest.constraints.IsEqualTo
import matlab.unittest.constraints.AbsoluteTolerance
import matlab.unittest.constraints.PublicPropertyComparator

this.verifyThat(actSolution, IsEqualTo(expSolution,...
    'Within', AbsoluteTolerance(10e-7),...
    'Using', PublicPropertyComparator.supportingAllValues));
end

author: translationalneuromodeling

repo: tapas