text:
function [x,f,fulltraining_instances,exitflag,output] = minFunc(funObj,x0,options,...
    params, allSNum, labels, cat_size, dictionary_length,freq_train, func, func_prime,We2,...
    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);

if isfield(options, 'logfile')
    logfile = options.logfile;
else
    logfile = [];
end

% 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] = feval(funObj, x, varargin{:});
else
    [f,g,H] = feval(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

if logfile
    fid = fopen(logfile, 'a');
    if (fid > 0)
        fprintf(fid, '-- %10s %10s %15s %15s %15s\n','Iteration','FunEvals','Step Length','Function Val','Opt Cond');
        fclose(fid);
    end
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
index=0;
% 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] = feval(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

    if logfile
        fid = fopen(logfile, 'a');
        if (fid > 0)
            fprintf(fid, '-- %10d %10d %15.5e %15.5e %15.5e\n',i,funEvals*funEvalMultiplier,t,f,sum(abs(g)));
            fclose(fid);
        end
    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;

    %% Fetch the features
    [W1, W2, W3, W4, b1, b2, b3, Wcat,bcat, We] = getW(1, x, params.embedding_size, cat_size, dictionary_length);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    disp('Get features by forward propagating and finding structure for all train and test sentences...')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % in this setting, we take the top node's vector and the average of all vectors in the tree as a concatenated feature vector
    fulltraining_instances = getFeatures(allSNum,0,...
        We,We2,W1,W2,W3,W4,b1,b2,b3,Wcat,bcat,params.alpha_cat,params.embedding_size, ...
        labels, freq_train, func, func_prime, params.trainModel);

    % 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: jacoxu


text:
function [M, P, B, A, PS] = matrix(disc, dim, domain)
%MATRIX    Convert operator to matrix using COLLOC discretization.
%   MATRIX(DISC) uses the parameters in DISC to discretize DISC.source as a
%   matrix using COLLOC.
%
%   MATRIX(DISC, DIM, DOMAIN) overrides the native 'dimension' and 'domain'
%   properties in DISC.
%
%   [PA, P, B, A, PS] = MATRIX(...) returns the projection matrix P, the
%   boundary matrix B, the unprojected cell array of square discretizations A,
%   and the projected conversion matrices PS.
%
% See also INSTANTIATE.

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

% Parse inputs
if ( nargin > 1 )
    disc.dimension = dim;
end
if ( nargin > 2 )
    disc.domain = domain;
end

% Check subinterval compatibility of domain and dimension.
if ( (length(disc.domain) - 1) ~= length(disc.dimension) )
    if ( length(disc.dimension == 1 ) )
        dim = disc.dimension;
        dom = disc.domain;
        disc.dimension = repmat(dim,1,length(dom)-1);
    else
        error('CHEBFUN:OPDISCRETIZATION:matrix:subIntDim', ...
            'Must specify one dimension value for each subinterval.')
    end
end

if ( nargout > 1 && ~isa(disc.source, 'linop') )
    error('CHEBFUN:OPDISCRETIZATION:matrix:matrix', ...
        'MATRIX() of a %s can only return one output.', class(disc.source))
end

if ( any(isinf(disc.domain)) )
    error('CHEBFUN:OPDISCRETIZATION:matrix:isinf', ...
        'Discretization on unbounded domains is not supported.');
end

% Construct a square representation of each block individually and
% store in a cell array.
[A, S] = instantiate(disc);

% We want output on different format depending on whether the source L is a
% LINOP or something else (typically a standard CHEBMATRIX):
if ( isa(disc.source, 'linop') )

    % Project rows down, and record the projection matrix as well.
    [PA, P, PS] = reduce(disc, A, S);

    % Get constraints:
    B = getConstraints(disc);

    % This should restore squareness to the final matrix.
    M = [ B ; PA ];

else

    % Everything should be of the same dimension.
    if ( iscell(A) )
        M = cell2mat(A);
    else
        M = A;
    end

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
clear;
load('../general/ccnf_patches_0.25_general.mat', 'centers', 'visiIndex', 'normalisationOptions');

mirrorInds = [1,17;2,16;3,15;4,14;5,13;6,12;7,11;8,10;18,27;19,26;20,25;21,24;22,23;...
  32,36;33,35;37,46;38,45;39,44;40,43;41,48;42,47;49,55;50,54;51,53;60,56;59,57;...
  61,65;62,64;68,66];

% For mirroring
frontalView = 1;
profileViewInds = [2,3,4];

% Grab all related experts and mirror them appropriatelly, just need to
% mirror the first layer

non_mirrored = mirrorInds(:,1);
normalisationOptions = rmfield(normalisationOptions, 'ccnf_ratio');
normalisationOptions.dccnf = true;

n_landmarks = size(visiIndex, 2);
n_views = size(visiIndex, 1);

patch_experts.correlations = zeros(n_views, n_landmarks);
patch_experts.rms_errors = zeros(n_views, n_landmarks);
patch_experts.types = {'reg'};
patch_experts.patch_experts = cell(n_views, n_landmarks);

scales = {'0.25', '0.35', '0.50', '1.00'};

root = 'D:/deep_experts/2017-02-02/rmses/';

for s=scales

    for c=1:n_views

        if(c == frontalView || sum(profileViewInds==c)> 0)

            for i=1:n_landmarks

                if(visiIndex(c,i))
                    mirror = false;
                    % Find the relevant file
                    if(c == frontalView)
                        rel_file = sprintf([root, 'MultiGeneral_arch4general_%s_frontal_%d_512.mat'], s{1}, i);
                    else
                        rel_file = sprintf([root, 'MultiGeneral_arch4general_%s_profile%d_%d_512.mat'], s{1}, c-1, i);
                    end
                    if(exist(rel_file, 'file'))
                       load(rel_file);
                    else
                       rel_id = mirrorInds(mirrorInds(:,2)==i,1);
                       if(isempty(rel_id))
                           rel_id = mirrorInds(mirrorInds(:,1)==i,2);
                       end
                       if(~visiIndex(c, rel_id))
                           continue;
                       end
                       if(c == frontalView)
                           rel_file = sprintf([root, 'MultiGeneral_arch4general_%s_frontal_%d_512.mat'], s{1}, rel_id);

                           mirror = true;
                           load(rel_file);
                       end
                    end
                    patch_experts.correlations(c, i) = correlation_2;
                    patch_experts.rms_errors(c, i) = rmse;

                    if(~mirror)
                        patch_experts.patch_experts{c, i} = weights;
                    else
                        flips = fliplr(reshape([1:121]', 11, 11));
                        weights_flipped = weights;
                        weights_flipped{1}(2:end,:) = weights{1}(flips+1,:);
                        patch_experts.patch_experts{c,i} = weights_flipped;
                    end
                end

            end
        else

            swap_id = find(centers(:,2) == -centers(c,2));

            corr_T = patch_experts.correlations(swap_id,:);
            % Appending a mirror view instead, based on the profile view
            corr_T  = swap(corr_T, mirrorInds(:,1), mirrorInds(:,2));
            patch_experts.correlations(c,:) = corr_T;

            rmsT = patch_experts.rms_errors(swap_id,:);
            rmsT = swap(rmsT, mirrorInds(:,1), mirrorInds(:,2));
            patch_experts.rms_errors(c,:) = rmsT;

            patchExpertMirror = patch_experts.patch_experts(swap_id,:);
            patchExpertMirrorT1 = patchExpertMirror(1,mirrorInds(:,1),:);
            patchExpertMirrorT2 = patchExpertMirror(1,mirrorInds(:,2),:);
            patchExpertMirror(1,mirrorInds(:,2),:) = patchExpertMirrorT1;
            patchExpertMirror(1,mirrorInds(:,1),:) = patchExpertMirrorT2;

            % To flip a patch expert it
            for p=1:size(patchExpertMirror,2)
                if(visiIndex(c, p))
                    weights = patchExpertMirror{p};
                    flips = fliplr(reshape([1:121]', 11, 11));
                    weights_flipped = weights;
                    weights_flipped{1}(2:end,:) = weights{1}(flips+1,:);
                    patch_experts.patch_experts{c,p} = weights_flipped;
                end
            end

        end
    end
    trainingScale = str2num(s{1});
    save(['cen_patches_', s{1} '_general.mat'], 'trainingScale', 'centers', 'visiIndex', 'patch_experts', 'normalisationOptions');

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function Tpix = conv_3d_T_from_phys_to_pix(Tphys, spc)
    Tpix = Tphys;
    for j = 1 : size(Tphys, 5)
    for i = 1 : size(Tphys, 4)
        Tpix(:,:,:, i, j) = Tpix(:,:,:,i, j) / spc(i);
    end
    end
end

repo: pTVreg


text:
% YALMIP
% Version 16-January-2020
% Help on http://yalmip.github.io

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%% Analyzing Neural Time Series Data
% Matlab code for Chapter 4 script C
% Mike X Cohen
%
% This code accompanies the book, titled "Analyzing Neural Time Series Data"
% (MIT Press). Using the code without following the book may lead to confusion,
% incorrect data analyses, and misinterpretations of results.
% Mike X Cohen assumes no responsibility for inappropriate or incorrect use of this code.

%%

% The function 'clear' removes all data from the matlab buffer
clear

% you can also clear specific variables.
leave_me_alone   = 10;
remove_me_please = 20;
clear remove_me_please
whos

remove_me_please % gives an error because the variable doesn't exist anymore!

%% basic importing text data

% The most basic way to import data is to copy and paste. This is the best
% option for small amounts of numeric data that you'll need to put into matlab
% only once. (Hint: use square brackets.)

% The simplest way to read in text data is if all data in the next are
% numbers (no text). Open a text editor and make a small matrix (say, 3x4).
% Next, type:
data = load('chapter04_datafile.txt');

% slightly more advanced:
[file_name,file_path]=uigetfile('*.txt'); % ui = user-interface
data = load([ file_path file_name ]);

% you can also read in data from excel files, but BE CAREFUL because this
% function can act in unexpected ways, e.g., by removing empty columns and
% rows without warning (this can be seen in comparing "numberdata" to "raw_data").
% Therefore, it might be best to use the "raw" data output.
[numberdata,textdata,raw_data] = xlsread('chapter04_excel_data.xls');

%% advanced importing text data

% Here we borrow from C language to flexibly read in mixed data. Let's say
% you have some poorly organized behavioral data files to read in, but at
% least you know what text strings to look for:

fid = fopen('chapter04_headache_data.txt','r');
% fid is a pointer to a location on the physical hard disk (similar to how
% we used variables as handles to axes when plotting). The 'r' means read
% (later we'll use 'w' for write).

% In this particular example, we will extract the trial number, subject
% choice, reaction time (RT), and accuracy for each trial. Fields are separated by tabs.

behavioral_data=[]; % initialize... we can't initialize the full matrix, because we don't know how big this will be.

% The following code will remain inside a loop, reading in and processing new
% lines of data, until we reach the end of the file.
datarow=1;

while ~feof(fid) % feof tests whether we're at the end of the file.

    dataline = fgetl(fid); % read a line ("file get line")

    dataline = regexp(dataline,'\t','split');
    % regexp can be used to cut data according to delimiters. Here we will
    % cut this string of characters into a cell array in which elements of
    % the array are separated by tabs.

    % here we use strcmpi to compare strings. The "i" means to ignore case.
    if ~any(strcmpi('trial',dataline))
        continue % continue means to skip to the next iteration of the loop.
    end

    trial_column    = find(strcmpi('trial',   dataline));
    choice_column   = find(strcmpi('choice',  dataline));
    rt_column       = find(strcmpi('rt',      dataline));
    accuracy_column = find(strcmpi('accuracy',dataline));

    behavioral_data(datarow,1) = str2double(dataline{trial_column+1});      % Note that we didn't initialize the size of the variable "behavioral_data" so matlab gives a warning.
    behavioral_data(datarow,2) = str2double(dataline{choice_column+1});     % If the variable is relatively small, it doesn't matter.
    behavioral_data(datarow,3) = str2double(dataline{rt_column+1});         % If the variable is large, however, it's best to initialize it to something really big, and then cut it down to size afterwards.
    behavioral_data(datarow,4) = str2double(dataline{accuracy_column+1});   % See chapter 4 in the book for further discussion of matrix initializations.

    datarow=datarow+1; % increment row
end

fclose(fid); % don't forget to close the file after you finish it!

%% initializing variables

num_rows = 10;
num_cols = 35;

% initialize with zeros (typical approach)
largematrix = zeros(num_rows,num_cols);

for rowi=1:num_rows
    for coli=1:num_cols
        % processing here...
    end % end row-loop
end % end column-loop

% note that you can increase the size of a matrix without initializing it
largematrix(num_rows+1,1) = 10;

% similarly, you can increase the dimensionality of a matrix
largematrix(1,round(num_cols/2),3) = 100;

% these last two options (adding elements and dimensions to an existing
% matrix) should be avoided whenever possible. They can create confusion and errors.

% you can also decrease matrix sizes/dimensions, ether by re-assignment:
largematrix = largematrix(:,:,1);

% or by setting parts of the matrix to be empty:
size(largematrix)
largematrix(:,end-4:end) = []; % this removes the last 5 columns
size(largematrix)

% Again, changing matrix sizes and dimensions should be avoided when
% possible, and done carefully when necessary.

%% basic saving data

% save as a .mat file (only matlab can read these files):
save('my_matlab_variables.mat','data','amsterdam','x'); % Question: Why does matlab crash on this line?

% The function 'dlmwrite' is useful if you have a matrix of numbers
% and want to write a text file of only numbers:
dlmwrite('data_written_from_matlab.txt',data,'\t');
% the final argument is the delimieter. This can be tab (\t), space ( ), comma (,), the letter X (X), etc.

%% advanced saving data

fid = fopen('data_output_SPSS_format.txt','w');

% we want the first row to be variable labels, then rows of mixed string-number data

% variable labels
variable_labels = {'Name';'trial';'choice';'rt';'accuracy'};

% let's add subject names
subject_names={'billy';'bob'};

for vari=1:length(variable_labels)
    fprintf(fid,'%s\t',variable_labels{vari});
    % the %s is for string; %g is for number.
end

% insert a new-line character
fprintf(fid,'\n');

for datarowi=1:size(behavioral_data,1)

    % print subject name
    fprintf(fid,'%s\t',subject_names{datarowi});

    % now loop through columns (variables)
    for columni=1:size(behavioral_data,2)
        fprintf(fid,'%g\t',behavioral_data(datarowi,columni));
    end
    fprintf(fid,'\n'); % end-of-line

    % You could also do this in one line:
    % fprintf(fid,'%s\t%g\t%g\t%g\t%g\n',subject_names{datarowi},behavioral_data(datarowi,1),behavioral_data(datarowi,2),behavioral_data(datarowi,3),behavioral_data(datarowi,4));

    fprintf('Finished writing line %g of %g\n',datarowi,size(behavioral_data,1));
end

fclose(fid);

% Now you can easily import these data into SPSS or Excel.

% Note that if you omit the first argument to fprintf, it puts the output
% in the command instead of the text file, as in the final line of this for-loop.

%% end.

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


text:
function [ y2, m2, d2, f2 ] = ymdf_inc_republican ( y1, m1, d1, f1, days )

%*****************************************************************************80
%
%% YMDF_INC_REPUBLICAN increments a Republican YMDF date by DAYS days.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 April 2103
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1,  D1, real F1,
%    the YMDF date.
%
%    Input, real DAYS, the number of days to advance the date.
%
%    Output, integer Y2, M2, D2, real F2,
%    the incremented YMDF date.
%

%
%  Copy the parameters.
%
  y2 = y1;
  m2 = m1;
  d2 = d1;
  f2 = f1 + days;
%
%  Check the parameters.
%
  [ y2, m2, d2, f2, ierror ] = ymdf_check_republican ( y2, m2, d2, f2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_bug1270

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY

% script to verify and fix an error related to using a mask in
% ft_singleplotER
%
% ??? Undefined function or variable 'xvector'.
%
% Error in ==> ft_singleplotER at 477
%    maskdatavector = reshape(mean(datmask,1), [1 numel(xvector)]);

clear all

mask = repmat([0 0 0 1 1 1 1 1 0 0],2,1);

timelock1 = [];
timelock1.label = {'1' '2'};
timelock1.time  = 1:10;
timelock1.dimord = 'chan_time';
timelock1.avg = randn(2,10);
timelock1.mask = timelock1.avg.*mask;

cfg = [];
cfg.channel = '1';
cfg.maskparameter = 'mask';
ft_singleplotER(cfg, timelock1);

% changing xvector into xval does the trick

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [f]= spm_fx_dem_pursuit(x,v,P)
% returns the flow for visual pursuit demo
% FORMAT [f]= spm_fx_dem_pursuit(x,v,P)
%
% x    - hidden states:
%   x.o(1) - oculomotor angle
%   x.o(2) - oculomotor angle
%   x.x(1) - target location (visual) - extrinsic coordinates (Cartesian)
%   x.x(2) - target location (visual) - extrinsic coordinates (Cartesian)
%   x.a(:) - attractor (SHC) states
%
% v    - hidden causes
% P    - parameters
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fx_dem_pursuit.m 4322 2011-05-04 15:28:08Z karl $

% intisaise flow (to ensure fields are aligned)
%--------------------------------------------------------------------------
f    = x;

% motion of attractor states
%==========================================================================
f.a  = spm_lotka_volterra(x.a,v);

% motion of target states
%==========================================================================

% target location is determined by the attractor state softmax(x.a)
%--------------------------------------------------------------------------
t    = P*spm_softmax(x.a,1/2);
f.x  = (t - x.x)/2;

% motion of oculomotor angles (attracted to target)
%==========================================================================
t    = atan(x.x);
f.o  = (t - x.o);

author: spm


text:
% Test file for trigtech/fliplr.m

function pass = test_fliplr(pref)

% Get preferences.
if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

%%
% Conduct a few very straightforward tests.

f = testclass.make(@(x) cos(pi*x), [], pref);
pass(1) = isequal(f, fliplr(f));

f = testclass.make(@(x) [sin(pi*x) cos(pi*x)], [], pref);
g = testclass.make(@(x) [cos(pi*x) sin(pi*x)], [], pref);
pass(2) = isequal(fliplr(f), g);

end


text:
function [LL, LLS, ht] = igarch_likelihood(parameters, epsilon, fepsilon, p, q, errorType, igarchType, constant, backCast, T, estimFlag)
% Log likelihood for IGARCH(P,Q) estimation
%
% USAGE:
%   [LL, LLS, HT] = igarch_likelihood(PARAMETERS,EPSILON,FEPSILON,P,Q,ERRORTYPE,IGARCHTYPE,CONSTANT,BACKCAST,T,ESTIMFLAG)
%
% INPUTS:
%   PARAMETERS   - A vector of IGARCH process parameters
%                    [omega alpha beta [nu lambda]]
%   DATA         - Vector of mean zero residuals
%   EPSILON      - A column of mean zero data
%   FEPSILON     - Either abs(EPSILON) or EPSILON.^2, depending on IGARCHTYPE
%   P            - Positive, scalar integer representing the number of
%                    symmetric innovations
%   Q            - Non-negative, scalar integer representing the number
%                    of lags of conditional variance (0 for ARCH)
%   ERRORTYPE    - [OPTIONAL] The error distribution used, valid types are:
%                    'NORMAL'    - Gaussian Innovations [DEFAULT]
%                    'STUDENTST' - T distributed errors
%                    'GED'       - Generalized Error Distribution
%                    'SKEWT'     - Skewed T distribution
%   IGARCHTYPE   - [OPTIONAL] The type of variance process, either
%                    1 - Model evolves in absolute values
%                    2 - Model evolves in squares [DEFAULT]
%   CONSTANT     - [OPTIONAL] Logical value indicating whether model
%                    should include a constant.  Default is true (include).
%   BACKCAST     - The value used for variance recursion
%   T            - Length of data
%   ESTIMFLAG    - [OPTIONAL] Flag (0 or 1) to indicate if the function
%                    is being used in estimation.  If it is 1, then the parameters are
%                    transformed from unconstrained values to constrained by standard
%                    garch model constraints
%
% OUTPUTS:
%   LL             - Minus 1 times the log likelihood
%   LLS            - Time series of log likelihoods (Also multiplied by -1)
%   HT             - Time series of conditional variances
%
% COMMENTS:
%   See also IGARCH

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 7/12/2009

if nargin==11 && estimFlag
    %If for estimation, transform the parameters
    [parameters,nu,lambda]=igarch_itransform(parameters,p,q,errorType,constant);
else
    %Otherwise the parameters simply must be parsed
    if errorType==2 || errorType==3
        %Seperate nu from the remaning parameters
        nu=parameters(p+q+constant);
        parameters=parameters(1:constant+p+q-1);
    elseif errorType==4
        lambda=parameters(p+q+constant+1);
        nu=parameters(p+q+constant);
        parameters=parameters(1:constant+p+q-1);
    end
end

%Backcast length
m  =  max([p q]);
%Compute the conditional variances
ht=igarch_core(fepsilon,parameters,backCast,p,q,m,T,igarchType,constant);
%Indices for the relevant opservations
t  = (m + 1):T;
ht = ht(t);
%Compute the log likelihoods
switch errorType
    case 1
        [LL, LLS] = normloglik(epsilon,0,ht);
        LLS = -LLS;
        LL = -LL;
    case 2
        [LL, LLS] = stdtloglik(epsilon,0,ht,nu);
        LLS = -LLS;
        LL = -LL;
    case 3
        [LL, LLS] = gedloglik(epsilon,0,ht,nu);
        LLS = -LLS;
        LL = -LL;
    case 4
        [LL, LLS] = skewtloglik(epsilon,0,ht,nu,lambda);
        LLS = -LLS;
        LL = -LL;
end

author: bashtage

repo: mfe-toolbox


text:
classdef CylindricalLiftModel < AbstractLiftCoefficientModel
    %CylindricalLiftModel Summary of this class goes here
    %   Detailed explanation goes here

    properties
        cylinderLength(1,1) double = 1; %m
        cylinderRadius(1,1) double = 1; %m

        liftCurves(1,1) LiftCoefficientCurves = LiftCoefficientCurves();
    end

    properties(Constant)
        enum = LiftCoefficientModelEnum.KSPCylinder;
    end

    methods
        function obj = CylindricalLiftModel(cylinderLength, cylinderRadius)
            obj.cylinderLength = cylinderLength;
            obj.cylinderRadius = cylinderRadius;

            obj.liftCurves = LiftCoefficientCurves();
        end

        function [ClS, liftUnitVectInertial] = getLiftCoeffAndDir(obj, ut, rVect, vVect, bodyInfo, mass, altitude, pressure, density, vVectECEF, attState)
            arguments
                obj(1,1) CylindricalLiftModel
                ut(1,1) double
                rVect(3,1) double
                vVect(3,1) double
                bodyInfo(1,1) KSPTOT_BodyInfo
                mass(1,1) double
                altitude(1,1) double
                pressure(1,1) double
                density(1,1) double
                vVectECEF(3,1) double
                attState(1,1) LaunchVehicleAttitudeState
            end

            [~,aoa,sideslip,totalAoA] = attState.getAeroAngles(ut, rVect, vVect, bodyInfo);

            vVectECEFMag = norm(vVectECEF);
            pressurePa = pressure*1000;
            speedSound = sqrt(1.4 * pressurePa / density); %m/s
            speedMS = vVectECEFMag*1000;                   %m/s
            thisMachNum = speedMS / speedSound;

            area = obj.getIncidentArea(totalAoA);

            ClS = obj.liftCurves.bodyLiftGiLift(sin(totalAoA))*obj.liftCurves.bodyLiftMachCurve(thisMachNum)*area;

            %get body centered coordinate systems
            bff = bodyInfo.getBodyFixedFrame();
            bci = bodyInfo.getBodyCenteredInertialFrame();

            %get the rotation matrix from body fixed to body centered
            %inertial
            R_ecef_to_global_inertial = bff.getRotMatToInertialAtTime(ut,[],[]);
            R_bci_to_global_inertial = bci.getRotMatToInertialAtTime(ut,[],[]);
            R_ecef_to_bci = R_bci_to_global_inertial' * R_ecef_to_global_inertial;

            %rotate the body X axis in the inertial frame to the body
            %fixed frame
            %TODO: Is attState.bodyX already in the body centered inertial
            %frame?  Or is it in the global inertial frame?
            R_bci_to_ecef = R_ecef_to_bci';
            bodyXInertial = attState.bodyX;
            bodyXECEF = R_bci_to_ecef*bodyXInertial;

            %get the lift vector in the body centered inertial frame
            v1 = crossARH(vVectECEF, bodyXECEF);
            liftUnitVectECEF = normVector(crossARH(v1, bodyXECEF));
            liftUnitVectInertial = R_ecef_to_bci * liftUnitVectECEF;
        end

        function useTf = openEditDialog(obj, lvdData)
            out = AppDesignerGUIOutput({false});
            lvd_EditKSPCylinderLiftPropertiesGUI_App(obj, lvdData, out);
            useTf = out.output{1};
        end
    end

    methods(Access=private)
        function area = getIncidentArea(obj, totalAoA)
            flatPlateArea = obj.cylinderLength*(2*obj.cylinderRadius);
            circleArea = pi*obj.cylinderRadius^2;

            area = cos(totalAoA)*circleArea + sin(totalAoA)*flatPlateArea;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [dLdp,iCpY,st] = spm_mci_grad_curve (assign,w,v,M,U,Y,fxtype)
% Compute gradient and curvature for MFX model
% FORMAT [dLdp,iCpY,st] = spm_mci_grad_curve (assign,w,v,M,U,Y,fxtype)
%
% assign    fields specify which are random/fixed effects
% w         random effects vector
% v         fixed effects vector
% M,U,Y     structure,inputs,data
% fxtype    'random' or 'fixed'
%
% dLdp      gradient
% iCpY      curvature (Fisher information)
% st        -1 for integration problem
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mci_grad_curve.m 6697 2016-01-27 14:57:28Z spm $

% Extract init and flow params from rfx or ffx vectors
[p_init,p_flow] = spm_mci_init_flow (assign,w,v,M);

Np=0;
st=0;

% Gradient and curvature of log likelihood
if strcmp(assign.init_par,fxtype)
    [G,sy_init,st] = spm_mci_sens_init (p_init,p_flow,M,U);
    sinit=1;
    Np=Np+size(sy_init,3);
else
    sinit=0;
end

if strcmp(assign.flow_par,fxtype)
    M.x0=p_init; % Initial conditions
    [G,sy_flow,st] = spm_mci_sens (p_flow,M,U);
    sflow=1;
    Np=Np+size(sy_flow,3);
else
    sflow=0;
end

if strcmp(assign.out_par,fxtype)
    M.x0=p_init; % Initial conditions
    % Compute sensitivity to output params
    [G,x,st] = spm_mci_fwd (p_flow,M,U);
    for n=1:M.N,
        [yout,dydx,dydoutp] = feval (M.g,x(n,:)',U(:,n),p_flow,M);
        sy_out(n,:,:)=dydoutp;
    end
    sout=1;
    Np=Np+size(sy_out,3);
else
    sout=0;
end

% Read data points and time indices
try, ind=Y.ind; catch, ind=1:M.N; end
Nt=length(ind);
y=Y.y;

if st==-1, disp('Integration Problem !'); return; end

% Prediction errors
g=G(ind,:);
e=Y.y-g;

% Compute gradient and precision
dLdp=zeros(1,Np);
iCpY=zeros(Np,Np);
for t=1:Nt,
    n=ind(t);
    sn=[];
    if sinit
        sn=squeeze(sy_init(n,:,:));
        if M.l==1, sn=sn'; end
    end
    if sflow
        sy=squeeze(sy_flow(n,:,:));
        if M.l==1, sy=sy'; end
        sn=[sn,sy];
    end
    if sout
        sy=squeeze(sy_out(n,:,:));
        if M.l==1, sy=sy'; end
        sn=[sn,sy];
    end
    dLdp=dLdp+e(t,:)*M.iCe*sn;
    iCpY=iCpY+sn'*M.iCe*sn;
end


text:
function [pval,rho] = SStat_mass_PHTest(X,d,t,Y,stats_PH,maskvtx)
% [pval,rho] = SStat_mass_PHTest(X,d,t,Y,stats_PH,maskvtx)
%
% Vertex/voxel-wise Schoenfeld residuals test for the proportional hazards
% assumption.
%
% Input
% X: Design Matrix with the time-independent covariates. (mxp, m # of
% subjects, p # of covariates).
% d: Logical vector (mx1) indicating censorship status (1 if the subject got
% the failure event or 0 otherwise).
% t: Vector (mx1) whose entries are the survival and censored times (ordered
% according to X).
% Y: Data matrix (mxnv, nv #vertices) whos colums represent the ordered
% data vector (according to X) for each voxel/vertex.
% stats_PH: Structure array containing statistics obtained from SStat_mass_CoxPH.
% maskvtx: Mask's vertices (1-based). Default [] (all vertices included).
%
% Output
% pval: Matrix of p-values (pxnv) indicating the probability of the PH
% assumption for each covariate and at each location.
% rho: Matrix of correlations (pxnv) between Schoenfeld residuals and ranked
% failure times at each location.
%
% $Revision: 1.1 $  $Date: 2015/01/06 17:03:59 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2015/01/06 17:03:59 $
%    $Revision: 1.1 $
% References: Kleinbaum, D.G., Klein, M., 2005. Survival analysis. A self-
% learning approach, second edition. New York: Springer..
%
%Validation of inputs
if nargin < 5
    error('Too few inputs');
elseif nargin < 6
    maskvtx = [];
end;
nv0 = size(Y,2);
if (length(stats_PH)~=nv0)
    error(['The structure array stats_PH must have the same length as the'...
        ' number of colums of Y.']);
end;
if isempty(maskvtx)
    maskvtx = 1:nv0;
end;
Y = Y(:,maskvtx);
nv = size(Y,2);
[m,p] = size(X);
if (length(d)~=m) || (length(t)~=m) || (size(Y,1)~=m)
    error(['The design matrix X, censorship status vector d, time'...
        ' vector t and data matrix Y must all have the same number of rows.']);
end;
%Sort the data by time. If there is a tie between a failure time and a
%censored time then the failure time goes first.
st_ix = find(d==1);
t1 = t(st_ix);
[t1,t1_ix] = sort(t1);
X1 = X(st_ix(t1_ix),:);
Y1 = Y(st_ix(t1_ix),:);
cs_ix = find(d==0);
if ~isempty(cs_ix)
    t2 = t(cs_ix);
    [t2,t2_ix] = sort(t2);
    X2 = X(cs_ix(t2_ix),:);
    Y2 = Y(cs_ix(t2_ix),:);
    count1 = 1; count2 = 1; i = 0;
    while (count1 <= length(t1)) && (count2 <= length(t2))
        i = i + 1;
        if t1(count1) <= t2(count2)
            X(i,:) = X1(count1,:);
            Y(i,:) = Y1(count1,:);
            d(i) = 1;
            t(i) = t1(count1);
            count1 = count1 + 1;
        else
            X(i,:) = X2(count2,:);
            Y(i,:) = Y2(count2,:);
            d(i) = 0;
            t(i) = t2(count2);
            count2 = count2 + 1;
        end;
    end;
    if (count1 > length(t1))
        X(i+1:end,:) = X2(count2:end,:);
        Y(i+1:end,:) = Y2(count2:end,:);
        d(i+1:end) = 0;
        t(i+1:end) = t2(count2:end);
    else
        X(i+1:end,:) = X1(count1:end,:);
        Y(i+1:end,:) = Y1(count1:end,:);
        d(i+1:end) = 1;
        t(i+1:end) = t1(count1:end);
    end;
else
    X = X1;
    Y = Y1;
    t = t1;
end;
%indices of unique failure times in ft_ix (last index when ties happen)
st_ix = find(d==1);
[ft,ft_ix] = unique(t(st_ix),'last');
ft_ix = st_ix(ft_ix);
nft = length(ft);
%handling ties in failure times by substracting a very small random number
rand('state',sum(100*clock));
for j=1:nft
    i = 1;
    while (ft_ix(j)-i>0) && (ft(j)==t(ft_ix(j)-i))
        i = i + 1;
    end;
    nties = i-1;
    tt = t(ft_ix(j)-nties:ft_ix(j)) - 1e-5*(1 + rand(nties+1,1));
    tX = X(ft_ix(j)-nties:ft_ix(j),:);
    tY = Y(ft_ix(j)-nties:ft_ix(j),:);
    [stt,stt_ix] = sort(tt);
    t(ft_ix(j)-nties:ft_ix(j)) = stt;
    X(ft_ix(j)-nties:ft_ix(j),:) = tX(stt_ix,:);
    Y(ft_ix(j)-nties:ft_ix(j),:) = tY(stt_ix,:);
end;
nft = length(st_ix);
p = p+1;
pval = ones(p,nv);
rho = zeros(p,nv);
r = zeros(nft,p);
display('Checking the proportional hazards assumption at each location ...');
for j=1:nv
    if ~isempty(stats_PH(j).Bhat)
        if length(stats_PH(j).Bhat)==p
            XD =  [X Y(:,j)];
            for i=1:nft
                lprv = exp(XD(st_ix(i):end,:)*stats_PH(j).Bhat);
                r(i,:) = XD(st_ix(i),:) - (lprv'*XD(st_ix(i):end,:))./sum(lprv);
            end;
            [rho(:,j),pval(:,j)] = corr(r,[1:nft]');
        else
            error([' The number of elements of stats(' num2str(j) ').Bhat'...
                ' is different from ' num2str(p)]);
        end;
    end;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function display(s)
%DISPLAY - spectralobject display overloaded operator

% VERSION: 1.0 of spectralobject
% AUTHOR: Celso Reyes
% LASTUPDATE: 1/30/2007

if isequal(get(0,'FormatSpacing'),'compact')
    disp([inputname(1) ' =']);
    disp(s);
else
    disp(' ');
    disp([inputname(1) ' =']);
    disp(' ');
    disp(s);
end

author: geoscience-community-codes

repo: GISMO


text:
function [ n_data, mu, sigma, a, x, fx ] = ...
  truncated_normal_a_pdf_values ( n_data )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_A_PDF_VALUES: values of the lower Truncated Normal AB PDF.
%
%  Discussion:
%
%    The Normal distribution, with mean Mu and standard deviation Sigma,
%    is truncated to the interval [A,+oo).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real MU, the mean of the distribution.
%
%    Output, real SIGMA, the standard deviation of the distribution.
%
%    Output, real A, the lower truncation limit.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 11;

  a_vec = [ ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0 ];

  fx_vec = [ ...
     0.01507373507401876, ...
     0.01551417047139894, ...
     0.01586560931024694, ...
     0.01612150073158793, ...
     0.01627701240029317, ...
     0.01632918226724295, ...
     0.01627701240029317, ...
     0.01612150073158793, ...
     0.01586560931024694, ...
     0.01551417047139894, ...
     0.01507373507401876 ];

  mu_vec = [ ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0 ];

  sigma_vec = [ ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0  ];

  x_vec = [ ...
     90.0, ...
     92.0, ...
     94.0, ...
     96.0, ...
     98.0, ...
    100.0, ...
    102.0, ...
    104.0, ...
    106.0, ...
    108.0, ...
    110.0 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    a = 0.0;
    mu = 0.0;
    sigma = 0.0;
    x = 0.0;
    fx = 0.0;
  else
    a = a_vec(n_data);
    mu = mu_vec(n_data);
    sigma = sigma_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function [logCount,t] = mandelbrotViewerProcessElement( x0, y0, escapeRadius2, maxIterations )
% Evaluate the Mandelbrot function for a single element

%   Copyright 2010-2011 The Mathworks, Inc.

t = 1;
z0 = complex( x0, y0 );
z = z0;
count = 0;
while count <= maxIterations && (z*conj(z) <= escapeRadius2)
    z = z*z + z0;
    count = count + 1;
end
magZ2 = max(real(z).^2 + imag(z).^2,escapeRadius2);
logCount = log( count + 1 - log( log( magZ2 ) / 2 ) / log(2) );

repo: Aero-matlab


text:
% save_intermediate - example hook function for each epoch
% Copyright (C) 2011 KyungHyun Cho, Tapani Raiko, Alexander Ilin
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%

function err = save_intermediate (R, args)

% args is a cell array
savefile = args{1};

save(savefile, 'R');

% -1 signals the end of training
err = 0;

author: kyunghyuncho

repo: deepmat


text: path1 = ''

author: xjqicuhk

repo: SIMS

save_path: github-repos/MATLAB/xjqicuhk-SIMS


text:
function w = ncc_compute_weights ( n )

%*****************************************************************************80
%
%% NCC_COMPUTE_WEIGHTS: weights of a Newton-Cotes Closed quadrature rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Output, real W(N), the weights.
%
  x_min = -1.0;
  x_max =  1.0;

  if ( n == 1 )

    w(1) = x_max - x_min;

  else

    x = linspace ( x_min, x_max, n );

    w = nc_compute ( n, x_min, x_max, x );

  end

  return
end


text:
function dawson_test ( )

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

  fprintf ( 1, '\n' );
  fprintf ( 1, 'DAWSON_TEST:\n' );
  fprintf ( 1, '  Test DAWSON_VALUES, R4_DAWSON, R8_DAWSON.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X      DAWSON(X)\n' );
  fprintf ( 1, '                 R4_DAWSON(X)         Diff\n' );
  fprintf ( 1, '                 R8_DAWSON(X)         Diff\n' );

  n_data = 0;

  while ( 1 )

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

    if ( n_data == 0 )
      break
    end

    fx2 = r4_dawson ( single ( x ) );
    fx3 = r8_dawson ( x );

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

  end

  rmpath ( '../test_values' );

  return
end


text:
function varargout = segGet(seg, property, varargin);
% Get properties of a segmentation.
%
% [vals] = segGet(seg, property, [optional arguments]);
%
% Some properties include:
%   'classPath':    path to the .class file describing the white/gray
%                   classification.
%   'grayPath':     path to the .gray graph.
%   'curMeshNum':   # of the currently-selected mesh (0=no meshes loaded).
%   'mesh',[M]:     the M-th loaded mesh <defaults to current mesh if M
%                   omitted.>
%
%   'nodes':        gray nodes. The rows of the 8xN nodes matrix represent:
%                   (1-3) s/i, a/p, l/r coords (in I|P|R space -- see
%                   mrFormatDescription) of each gray node.
%                   (4)
%                   (5)
%                   (6) gray level/layer, counting from white matter up.
%                   (7)
%                   (8)
%                   The number N of nodes will represent all the gray
%                   nodes for this segmentation, regardless of the
%                   current mapping.
%   'edges':        gray edges.
%   [nodes edges] = segGet('gray') returns both gray nodes and edges
%                   (preventing redundant loading of the gray graph).
%
%   'meshCoords',[M],[mmPerPix]:
%                   3xN coordinates (in IPR space) of each node
%                   represented by the M-th mesh, using the current
%                   mapping. The current mapping is determined by
%                   mrmPreferences. The number N of columns will
%                   correspond to the columns of mesh{M}.initVertices,
%                   as well as the mesh colors.
%                   <if M omitted, uses selected mesh; if mmPerPix
%                   omitted, tries to read from segementation's anat file,
%                   or else assumes 1x1x1.>
%
%	'nearestNode',[coord],[tolerance=5mm]:
%					Return the index of the nearest gray node to a given
%					3D coordinate. This simply finds the node with the
%					smallest euclidean distance from the input coordinate.
%					The input coord should be a 1x3 or 3x1 vector, and should
%					be in the same format as the gray coords. (This is in order
%					The optional tolerance argument specifies how large a
%					distance should be accepted for a nearest node; if no
%					node is less than this amount, the nearestNode will
%					return empty. By default, the tolerance is 5mm.
%
% ras 04/06
if nargin<2, help(mfilename); error('Not enough input args.'); end

switch lower(property)
    case {'classpath' 'classfile'}
        if isstr(seg.class),    varargout{1} = seg.class;
        else,                   varargout{1} = seg.class.filename;
        end

    case {'classification' 'classdata' 'voxels' 'wm' 'class'}
        if isstr(seg.class),    varargout{1} = readClassFile(seg.class);
        else,                   varargout{1} = seg.class;
        end

    case {'graypath' 'graygraph' 'grayfile'}
        if isstr(seg.gray),     varargout{1} = seg.gray;
        else,                   varargout{1} = seg.gray.path;
        end

    case {'curmeshnum' 'curmeshn' 'meshnum' 'selectedmeshnum'}
        varargout{1} = seg.settings.mesh;

    case {'mesh' 'curmesh' 'selectedmesh'}
        if length(varargin)==0, M = segGet(seg, 'curMeshNum');
        else, M = varargin{1};
        end

        if M==0 | isempty(seg.mesh),
            varargout{1} = [];
        else
            varargout{1} = seg.mesh{M};
        end

    case {'nodes' 'graynodes'}
        if isempty(seg.nodes),  varargout{1} = readGrayGraph(seg.gray);
        else,                   varargout{1} = seg.nodes;
        end

    case {'edges' 'grayedges'}
        if isempty(seg.edges),  [x varargout{1}] = readGrayGraph(seg.gray);
        else,                   varargout{1} = seg.edges;
        end

    case {'gray'}
        if isempty(seg.nodes) | isempty(seg.edges)
            [varargout{1} varargout{2}] = readGrayGraph(seg.gray);
        else
            varargout{1} = seg.nodes;
            varargout{2} = seg.edges;
        end

    case {'graycoords'}
        nodes = segGet(seg, 'nodes');
        varargout{1} = nodes([2 1 3],:);

    case {'meshcoords' 'coords'}
        varargin{3} = []; % fast way to init. args 1 & 2 if unspecified
        M = varargin{1};  mmPerPix = varargin{2};
        if isempty(M), M = segGet(seg, 'curMeshNum'); end
        if isempty(mmPerPix)
            try
                anat = mrLoad(seg.anatFile);
                mmPerPix = anat.voxelSize(1:3);
            catch
                disp(['Couldn''t read segmentation anatomy. Guessing ' ...
                      'mmPerPix as [1 1 1]...'])
                 mmPerPix = [1 1 1];
            end
        end

        % get mapping from mesh vertices to gray nodes
        [nodes edges] = segGet(seg, 'gray');
        v2g = mrmMapVerticesToGray(seg.mesh{M}.initVertices, nodes, ...
                                   mmPerPix, edges);

        % grab appropriate coords from nodes
        varargout{1} = nodes(1:3,v2g);

	case {'nearestnode'}
        varargin{3} = []; % fast way to init. args 1 & 2 if unspecified
		pt = varargin{1}(:); tolerance = varargin{2};
		if isempty(pt), error('Need an input coordinate.'); end
		if isempty(tolerance), tolerance = 5; end

		% compute Euclidean distance of each node from the coord
		C = segGet(seg, 'GrayCoords');
		dist = sqrt( [C(1,:) - pt(1)] .^ 2 + ...
					 [C(2,:) - pt(2)] .^ 2 + ...
					 [C(3,:) - pt(3)] .^ 2 );
		if min(dist) > tolerance
			warning(sprintf('[%s]: no node found within tolerance [%i mm]', ...
							 mfilename, tolerance));
			varargout{1} = [];
		else
			I = find(dist==min(dist));
			varargout{1} = I(1);
		end

    otherwise, error('Unknown property.')

end

return


text:
% OtSkimCubeNew Create a new skim cube.
%
% Usage:
%    hscube = OtSkimCubeNew(<path to a work directory>)
% or:
%    hscube = OtSkimCubeNew()
%    to try to use a temporary directory (%TMP%, %TEMP% or /tmp)
%
% To release a skim cube handle, call OtSkimCubeFree().

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-04-19 Created (Johan Meijdam, VORtech)

function hscube = OtSkimCubeNew(workdir)
   if nargin == 0
      workdir = getenv('TMP');
      if isempty(workdir)
         workdir = getenv('TEMP');
         if isempty(workdir)
            workdir = [filesep 'tmp'];
         end
      end
   end
   hscube = Ot('skimcube.new', workdir);
end

author: Sable

repo: mcbench-benchmarks


text:
function [bboxes,idx] = doNMS(bboxes, thresh)

if isempty(bboxes)
    return
end
tmp = bboxes(:,1);
idx = 1;
for i = 2:size(bboxes,2)
    if max(getIOUFloat(tmp',bboxes(:,i)')) < thresh
        tmp = [tmp bboxes(:,i)];
        idx = [idx i];
    end
end
bboxes = tmp;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function geometry_test050 ( )

%*****************************************************************************80
%
%% TEST050 tests PLANE_EXP_NORMAL_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST050\n' );
  fprintf ( 1, '  PLANE_EXP_NORMAL_3D finds the normal \n' );
  fprintf ( 1, '    to a plane.\n' );

  p1(1:dim_num) = [ -10.56, -10.56, 78.09 ];
  p2(1:dim_num) = [  44.66, -65.77,  0.00 ];
  p3(1:dim_num) = [  44.66,  44.66,  0.00 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Three points on the plane:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  P1: %10f  %10f  %10f\n', p1(1:dim_num) );
  fprintf ( 1, '  P2: %10f  %10f  %10f\n', p2(1:dim_num) );
  fprintf ( 1, '  P3: %10f  %10f  %10f\n', p3(1:dim_num) );

  normal = plane_exp_normal_3d ( p1, p2, p3 );

  r8vec_print ( dim_num, normal, '  The normal vector:' );

  return
end


text:
function checksum=NMEAChecksum(NMEAString)
%%NMEACHECKSUM Compute the checksum of a National Maritime Electronics
%              Association (NMEA) 0183 sentence. This is the message
%              standard used for Automatic Identification System (AIS)
%              messages sent by transponders on ships. This function is
%              useful for validating checksums to determine whether to
%              discard a message as corrupt.
%
%INPUTS: NMEAString A NMEA 0183 string as a sequence of characters. The
%                   string should start with a ! or $ and everything until
%                   a * goes into the checksum. In such a message, the
%                   checksum is placed after the *. Additional data might
%                   come after that.
%
%OUTPUTS: checksum A two-character checksum associated with the NMEA
%                  string. Everything after the * is ignored when
%                  this function computes the checksum. If the message
%                  does not have characters surrounded by a ! or $ and
%                  a *, then an empty matrix is returned. The function does
%                  no further checks to make sure that the message is in
%                  the correct format.
%
%The NMEA 0183 standard is expensive and can be purchased from
%http://www.nmea.org/content/nmea_standards/nmea_0183_v_410.asp
%However, a number of web pages say that the checksum is just the
%hexadecimal representation of xoring the ASCII codes of the characters
%between the ! or $ and the * in the message.
%
%As an example, calling the function as
%NMEAChecksum('!AIVDM,1,1,,B,177KQJ5000G?tO`K>RA1wUbN0TKH,0*5C')
%where the message is an AIS message would return '5C", indicating that the
%checksum embedded in the message is correct. Alternatively, if one were
%trying to create a checksum to be appended to the message, one would call
%NMEAChecksum('!AIVDM,1,1,,B,177KQJ5000G?tO`K>RA1wUbN0TKH,0*')
%and would get the same result as everything after the * is ignored.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%If it has an obviously wrong format.
if(isempty(NMEAString)||~ischar(NMEAString)||(~isrow(NMEAString)&&~iscolumn(NMEAString))||(NMEAString(1)~='$'&&NMEAString(1)~='!')||length(NMEAString)<3)
    checksum=[];
    return;
end

%Find the first * in the string.
maxIdx=2;
while(NMEAString(maxIdx)~='*')
    maxIdx=maxIdx+1;

    if(maxIdx>length(NMEAString))
        checksum=[];
        return;
    end
end
maxIdx=maxIdx-1;

%The bitxor function will only work with raw UTF8/ASCII codes.
NMEAString=unicode2native(NMEAString);

checksum=0;
for curIdx=2:maxIdx
    checksum=bitxor(checksum,NMEAString(curIdx));
end
checksum=dec2hex(checksum,2);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [elem,idx,area,bdFlag] = fixorder(node,elem,bdFlag)
%% FIXORDER set all triangles counter-clockwise
%
%   elem = FIXORDER(node,elem) computes signed area (volume) of all
%   triangles (tetrahedron) in the triangulation and switch the vertices
%   such that all signed area (volume) are positive.
%
%   [elem,idx,area] = FIXORDER(node,elem) also outputs the index set of
%   elements whose area (volume) are negative and the absolute value of area.
%
%   [elem,idx,area,bdFlag] = FIXORDER(node,elem,bdFlag) changes the bdFlag
%   for boundary conditions.
%
%   fixorder is recommend to use if the mesh is obtained by delaunay
%   function in matlab.
%
% See also fixorientation, fixorder3
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if (nargin==2)
    bdFlag = [];
end
% compute signed area of each triangle
[area,elemSign] = simplexvolume(node,elem);
% find triangles with negative area and switch the vertices
idx = find(elemSign==-1);
elem(idx,[2 3]) = elem(idx,[3 2]);
if exist('bdFlag','var') && ~isempty(bdFlag)
    bdFlag(idx,[2 3]) = bdFlag(idx,[3 2]);
end

repo: ifem


text:
function theta = get_theta_bounds(e1,e2,F,dim)

[theta1,~] = get_theta(e1,dim);
[theta2,l2] = get_theta(e2,dim);

%% trasfer l2 to i1
l2(:,1) = F_transfer_l(F,l2(:,1),e2);
l2(:,2) = F_transfer_l(F,l2(:,2),e2);

theta2(1) = l_get_angle(l2(:,1));
theta2(2) = l_get_angle(l2(:,2));

if theta2(1)>theta2(2)
    theta2 = flipud(theta2);
end

theta(1) = max(theta1(1),theta2(1));
theta(2) = min(theta1(2),theta2(2));


text:
function [tBreak,a,b,c,m,p]=essential_distMinAnglePair_discontinuityDistance(Q21)
a=Q21(1,1)+Q21(2,2);
b=Q21(1,2)-Q21(2,1);
c=Q21(3,3);

m=norm([a;b]);
p=sign(a)*acos(clip(b/m));

%tBreak=modAngle(3/2*pi-p);
tBreak=-0.5*pi-p;

function v=clip(v)
v=min(1,max(-1,v));

author: NicolasBoumal

repo: manopt