text:
function [D, S] = spm_eeg_reref_eeg(S)
% Rereference EEG data to new reference channel(s)
% FORMAT [D, S] = spm_eeg_reref_eeg(S)
%
% S                    - input structure (optional)
% (optional) fields of S:
%   S.D                - MEEG object or filename of M/EEG mat-file
%   S.refchan          - New reference channel indices or labels
%                        ('average' can be used as shortcut)
%
% Output:
% D                    - MEEG object (also written on disk)
% S                    - record of parameters, including montage
%__________________________________________________________________________
%
% spm_eeg_reref_eeg re-references any EEG data within an MEEG dataset, by
% calling spm_eeg_montage with appropriate montage, excluding bad channels
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Rik Henson (updated so that bad channels included in tra but not used)
% $Id: spm_eeg_reref_eeg.m 6374 2015-03-12 10:15:18Z vladimir $

SVNrev = '$Rev: 6374 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','EEG reref'); spm('Pointer','Watch');

%-Get MEEG object
%--------------------------------------------------------------------------
try
    D = S.D;
catch
    [D, sts] = spm_select(1, 'mat', 'Select M/EEG mat file');
    if ~sts, D = []; return; end
    S.D = D;
end

D = spm_eeg_load(D);

% Get indices for just EEG channels and remove any bad channels
%--------------------------------------------------------------------------
eegchan  = D.indchantype('EEG');
goodchan = setdiff(eegchan, D.badchannels);

if isempty(goodchan)
    error('No EEG channels that are not marked as bad...!?')
end

%-Get reference channel indices
%--------------------------------------------------------------------------
if ~isfield(S, 'refchan')
      [selection, ok]= listdlg('ListString', D.chanlabels(goodchan), 'SelectionMode', 'multiple' ,'Name', 'Select reference channels' , 'ListSize', [400 300]);
    if ~ok
        return;
    end

    if length(selection) == length(goodchan)
        S.refchan = 'average';
    else
        S.refchan = D.chanlabels(goodchan(selection));
    end
end

if strcmpi(S.refchan,'average')
    refchan = goodchan;
elseif iscell(S.refchan)
    refchan = D.indchannel(S.refchan);
elseif isnumeric(S.refchan)
    refchan = S.refchan;
end

refind  = find(ismember(eegchan, refchan));
goodind = find(ismember(eegchan, goodchan));
badind  = find(ismember(eegchan, D.badchannels));

if length(refind) ~= length(refchan)
    error('Not all S.refchan are valid indices for (non-bad) EEG channels')
end

tra                 = eye(length(eegchan));
tra(goodind,refind) = tra(goodind,refind) - 1/length(refchan);
tra(badind,refind)  = tra(badind,refind)  - 1/length(refchan);

S1=[];
S1.D = D;
S1.montage.labelorg = D.chanlabels(eegchan);
S1.montage.labelnew = D.chanlabels(eegchan);
S1.montage.chantypeorg = lower(D.chantype(eegchan)');
S1.montage.chantypenew = lower(S1.montage.chantypeorg);
S1.montage.chanunitorg = D.units(eegchan)';
S1.montage.chanunitnew = S1.montage.chanunitorg;
S1.montage.tra = tra;
S1.keepothers = 1;
S1.updatehistory = 0;
D = spm_eeg_montage(S1);

%-Update history (not necessary; leave as call to spm_eeg_montage?) Save new evoked M/EEG dataset
%--------------------------------------------------------------------------
D = D.history('spm_eeg_reref_eeg', S);

save(D);

%-Cleanup
%--------------------------------------------------------------------------
spm_progress_bar('Clear');
spm('FigName','EEG re-reference: done'); spm('Pointer','Arrow');

author: spm


text:
function avw = ctf_mri2avw(mri)

% avw = ctf_mri2avw(mri)
%
% The purpose of this function is to convert a CTF .mri volume into
% an Analyze volume.  The returned avw struct can be saved using avw_write.
%
% This function depends on the avw* functions availabe in the mri_toolbox
% at http://eeg.sf.net/.
%
%      <>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> %
%      <                                                      > %
%      <                    DISCLAIMER:                       > %
%      <                                                      > %
%      < THIS PROGRAM IS INTENDED FOR RESEARCH PURPOSES ONLY. > %
%      < THIS PROGRAM IS IN NO WAY INTENDED FOR CLINICAL OR   > %
%      <                    OFFICIAL USE.                     > %
%      <                                                      > %
%      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<> %
%

% $Revision: 1.1 $ $Date: 2009-01-30 03:49:27 $

% Copyright (C) 2004  Darren L. Weber
%
% 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

% History:  04/2004, Darren.Weber_at_radiology.ucsf.edu
%                    - adapted from an appendex to CTF document
%                    MRIConverter.pdf, which is copied at the end of this
%                    function.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

avw = avw_hdr_make;

avw.fileprefix = strrep(mri.file,'.mri','');

ver = '[$Revision: 1.1 $]';
fprintf('\nCTF_MRI2AVW [v%s]\n',ver(12:16));  tic;

% CTF .mri files are always 256x256x256 voxels, 1mm^3
avw.hdr.dime.dim(2) = 256;
avw.hdr.dime.dim(3) = 256;
avw.hdr.dime.dim(4) = 256;
avw.hdr.dime.pixdim(2) = 1;
avw.hdr.dime.pixdim(3) = 1;
avw.hdr.dime.pixdim(4) = 1;

% mri.hdr.dataSize = 1 or 2 (bytes), 8 or 16 bits
if mri.hdr.dataSize == 1,
    avw.hdr.dime.bitpix = 8;
    avw.hdr.dime.datatype = 2;
else
    avw.hdr.dime.bitpix = 16;
    avw.hdr.dime.datatype = 4;
end

% This next step should always work correctly, given that avw structs
% are always axial unflipped in the matlab workspace; the axial
% unflipped orientation is
% (+X left,  +Y anterior,  +Z superior); the ctf orientation is
% (+X right, +Y posterior, +Z inferior), the complete opposite.
temp = flipdim(mri.img,1);
temp = flipdim(temp,2);
temp = flipdim(temp,3);
avw.img = temp;

t=toc; fprintf('...done (%5.2f sec).\n\n',t);

return

author: PatternRecognition

repo: OpenBMI


text:
function [xk,objfunvals,xkk] = solveLASSOProblem(A,b,alpha,varargin)
%%SOLVELASSOPROBLEM This function solves the Least Absolute Shrinkage and
%               Selection Operator (LASSO) problem. This optimization
%               problem is
%               minimize (1/2)*norm(A*x-b,2)^2+alpha*norm(x,1)
%               Thus, the LASSO problem is the same as the ridge regression
%               problem except the second term uses an l_1 norm instead of
%               an l_2 norm. This type of optimization plays a notable role
%               in compressive sensing algorithms. The optimization problem
%               is solved using the Fast-Iterative Shrinkage and
%               Thresholding Algorithm (FISTA) for sparse vectors.
%
%INPUTS:      A An mXn matrix. This can be complex.
%             b An mX1 vector. This can be complex.
%         alpha A real scalar value.
%      varargin Arguments that are specified as fieldname, value pairs.
%               Possible pairs are:
%               FIELD NAME         VALUE TYPE     DESCRIPTION
%               -----------------  -------------- -----------------------
%               'maxiter'          numeric scalar The maximum number of
%                                                 iterations (default 2000)
%               'tol'              numeric scalar Tolerence for termination
%                                                 (default 1e-6).
%               'tau'              numeric scalar Descent parameter used to
%                                                 to control the rate of
%                                                 minimization (default:
%                                                 svds(A,1).^2).
%               'backtracking'     bool           Performs backtracking if
%                                                 turned on.
%               'verbose'          bool           Enables verbose mode.
%               'stoppingcriteria' char           'cauchy' or 'gradient'.
%               'init'             vector         n x 1 vector for warm
%                                                 start.
%
%OUTPUTS: xk The nX1 vector solution to the optimization problem.
% objfunvals If verbose mode is activated, then this will hold the value of
%            the objective function at every iteration. in a numIterX1
%            vector. Otherwise, this is an empty matrix.
%        xkk An nXnumIter matrix holding the value of the estimate at each
%            step.
%
%The algorithm implemented is that of [1]. The algorithm typically needs a
%very large number of iterations to converge. However, the number of
%iterations is not particularly dependent on the scale of the problem,
%making the solver good for large problems.
%
%EXAMPLE:
%A trivial example is
% A=[10,  4,   -3, 6;
%     0,  -18, 24, -10;
%     5,  7,   3,  12];
% xTrue=[0;2;0;12];
% alpha=2;
% b=A*xTrue;
% xk=solveLASSOProblem(A,b,alpha)
%The xk found should have the required sparsity, but the l1 penalty term
%means that the xk found will not equal xTrue; Rather, it will be
%xk=[0;2.006;0;11.9922]. Increasing alpha forces the sparsity up. For
%example, using alpha=3200 changes the result to xk=[0;2.6581;0;0] and
%using alpha=4000 changes xk to all zeros. Making alpha close to 0 reduces
%the result to approach pinv(A)*b.
%
%REFERENCES:
%[1] A. Beck and M. Teboulle, "A fast iterative shrinkage-thresholding
%    algorithm for linear inverse problems," SIAM Journal on Imaging
%    Sciences, vol. 2, no. 1, pp. 183-202, 2009.
%
%February 2016 Hatim F. Alqadah, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    %% PROCESS INPUTS HERE
    if(nargin < 3)
        alpha = [];
    end
    if(nargin < 2)
        error('invalid number of arguments');
    end

    % check argument type and sizes
    assert(ismatrix(A),'A should be a matrix');
    assert(isvector(b),'b should be a vector');
    [m,n] = size(A);
    assert(size(b,1) == m,'invalid vector length on right hand side');

    % if matrix or data is complex we need to seperate out
    % the real and imaginary components
    is_complex = ~(isreal(A) && isreal(b));
    if(is_complex)
        % reformulate into a purely real problem
        A = [real(A) -imag(A);imag(A) real(A)];
        b = [real(b);imag(b)];
    end

    % set default options and process user defined ones
    opts.maxiter = 2000;
    opts.tol = 1e-6;
    opts.stopcrit = 'cauchy';
    opts.backtracking = false;
    opts.verbose = false;
    opts.restart = false;
    opts.return_as_dense = true;
    opts.init = [];
    opts.tau = [];

    [opts] = addCelListToStruct(opts,varargin);
    clear varargin

    if(~isempty(opts.init))
        assert(length(opts.x0) == n,'given warm start has incorrect length');
        xk = x0;
        if(is_complex)
            xk = [real(xk);imag(xk)];
        end
    else
        if(is_complex)
            xk = sparse(2*n,1);
        else
            xk = sparse(n,1);
        end
    end

    %% PREPROCESSING

    % set initial quantities
    xp = xk;
    nn = 0;
    tk = 1;
    tp = 1;
    At = A';
    Atb = At*b;
    has_converged = false;
    keep_history = nargout > 2;
    objfunvals = [];
    xkk = [];

    if(isempty(alpha))
        alpha = .5*norm(Atb,'inf');
    end

    if(isempty(opts.tau))
        %Use iterative method to find the Lipschitz constant
        opts.tau = svds(A,1).^2;
    end

    %Setup reporting functions if quiet mode is disabled
    if(opts.verbose)
        fprintf('\n*******************************************\n');
        fprintf('* FISTA: SPARSE VECTOR RECONSTRUCTION \n');
        if(is_complex)
        	fprintf('* CMPLX     : %s\n','true');
        else
            fprintf('* CMPLX     : %s\n','false');
        end
        fprintf('* MAXITER   : %d\n',opts.maxiter);
        fprintf('* TOL       : %e\n',opts.tol);
        if(opts.restart)
            fprintf('* RESTART   : %s\n','true');
        else
            fprintf('* RESTART   : %s\n','false');
        end
        fprintf('* STOPCRIT  : %s\n',opts.stopcrit);
        fprintf('*******************************************\n');
        msgstr = ['Current Precision: ',num2str(inf)];
        h = waitbar(0,msgstr);
        precision = zeros(opts.maxiter,1);
        objfun = @(x) ObjFunc(x,A,b,alpha);
        objfunvals = zeros(opts.maxiter,1);
    end
    if(keep_history)
        xkk = zeros(length(xk),opts.maxiter);
    end

    %% MAINLOOP

    while(nn < opts.maxiter)
        if opts.verbose
            objfunvals(nn+1) = objfun(xk); %#ok<*AGROW>
        end
        if(keep_history)
            xkk(:,nn+1) = xk;
        end

        yk = xk + ((tp-1)/tk)*(xk - xp);
        grd = At*(A*yk) - Atb;
        xp = xk;
        if(opts.backtracking)
            % solve the local regularization problem with an appropriate
            % step size
            [xk,opts.tau] = LineSearch(A,b,grd,alpha,yk,opts.tau);
        else
            xk = sparse(shrink(yk - (1/opts.tau)*grd,alpha/opts.tau));
        end

        % check stopping criterion
        switch(lower(opts.stopcrit))
            case 'cauchy'
                prec = norm(xk(:) - yk(:))/max(norm(xk),1);
            case 'subgradient'
                % we need to estimate the subgradient at current pt
                % optimality occurs at zero
                prec = ComputeOptimalityResidue(A,At,xk,Atb,alpha);
        end
        if(prec < opts.tol)
            has_converged = true;
            break;
        end

        % keep going
        nn = nn + 1;
        tp = tk;
        tk = 0.5*(1+sqrt(1+4*tk*tk));

        % check if a restart is necessary
        if(opts.restart && ((yk - xk)'*(xk - xp) > 0))
            % forget all previous iterations
            % and reset momentum back to zero
            tp = 1;
            tk = 1;
        end
        if(opts.verbose)
            precision(nn) = prec;
            waitbar(nn/opts.maxiter,h,['Current precision: ',num2str(prec)]);
        end
    end
    if(~has_converged)
        warning('MATLAB:no_convergence','method failed to converge within given time steps');
    end

    %% Post-processing
    if(opts.verbose)
        waitbar(1,h);
        objfunvals = objfunvals(1:nn);
        precision = precision(1:nn);
        iter = 0:nn-1;
        figure;
        semilogy(iter,objfunvals,'--r','linewidth',1.2);
        grid on;
        title(['Objective Function for \alpha = ',num2str(alpha)]);
        xlabel('ITERATION');
        ylabel('FOBJ VAL');

        figure;
        semilogy(iter,precision,'--g','linewidth',1.2);
        grid on;
        title(['Precision for \alpha = ',num2str(alpha)]);
        xlabel('ITERATION');
        ylabel('PREC LVL');
        close(h);
    end
    if(opts.return_as_dense)
        xk = full(xk);
        if keep_history
            xkk = full(xkk(:,1:nn));
        end
    end
    if(is_complex)
        xk = complex(xk(1:n),xk(n+1:end));
        if(keep_history)
            xkk = complex(xkk(1:n,:),xkk(n+1:end,:));
        end
    end

end
%% UTILITY
% defines the standard shrinkage operator.
function z = shrink(x,mu)
    z = sign(x).*max(abs(x)-mu,0);
end

% Experimental stopping criteria particlarly for proximal gradient methods.
% It is not clear if this can help issues do non-strict convexivity of the
% data fidelity term
function prec = ComputeOptimalityResidue(A,At,xk,Atb,alpha)
    % first compute the gradient at the current point
    gradF = At*(A*xk) - Atb;

    % compute the unique portion of the sub-gradient
    subGrad = alpha*sign(xk);

    % now for the points that are zero minimize the
    % optimality residue
    inds = find(subGrad == 0);

    % compute solution to the constrained
    % sub optimization problem
    g = gradF(inds);
    y = zeros(size(g));

    for i=1:50
        yp = y;
        lambda = abs(y+g)/alpha;
        y = (y - .5*(y+g))./(1 - .5*lambda);
        if(norm(y-yp)/max([norm(y),1]) < 1e-5)
            break;
        end
    end
    subGrad(inds) = y;

    % evaluate the maximum sub-gradient component
    prec = norm(gradF + subGrad,'inf');
end

% returns the objective function value of the unconstrained l2-l1 problem
function val = ObjFunc(x,A,b,alpha)
   val = .5*norm(A*x - b).^2 + alpha*norm(x,1);
end

% A straight-forward line search method for first order methods
function [xk,fxk,L] = LineSearch(A,b,grd,alpha,yk,L)
    beta = 1.2;
    stopBacktrack = false;
    fyk = .5*norm(A*yk - b).^2 + alpha*norm(yk,1);
    while(~stopBacktrack)
        % predict the next value
        xk = sparse(shrink(yk - (1/L)*grd,alpha/L));
        fxk = .5*norm(A*xk - b).^2 + alpha*norm(xk,1);

        % find the Lipschitz constant that ensures the Taylor
        % expansion is valid
        q =  fyk + (xk-yk)'*grd + 0.5*L*norm(xk-yk).^2 + alpha*norm(xk,1);
        if(fxk <= q)
            stopBacktrack = true;
        else
            % increase the Lipschitz constant estimate
            L = L*beta;
        end
    end
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 [t,x_new,f_new,g_new,funEvals,H] = ArmijoBacktrack(...
    x,t,d,f,fr,g,gtd,c1,LS_interp,LS_multi,progTol,debug,doPlot,saveHessianComp,funObj,varargin)
% [t,x_new,f_new,g_new,funEvals,H] = ArmijoBacktrack(...
%    x,t,d,f,fr,g,gtd,c1,LS_interp,LS_multi,progTol,debug,doPlot,saveHessianComp,funObj,varargin)
%
% Backtracking linesearch to satisfy Armijo condition
%
% Inputs:
%   x: starting location
%   t: initial step size
%   d: descent direction
%   f: function value at starting location
%   fr: reference function value (usually funObj(x))
%   gtd: directional derivative at starting location
%   c1: sufficient decrease parameter
%   debug: display debugging information
%   LS_interp: type of interpolation
%   progTol: minimum allowable step length
%   doPlot: do a graphical display of interpolation
%   funObj: objective function
%   varargin: parameters of objective function
%
% Outputs:
%   t: step length
%   f_new: function value at x+t*d
%   g_new: gradient value at x+t*d
%   funEvals: number function evaluations performed by line search
%   H: Hessian at initial guess (only computed if requested)
%
% recet change: LS changed to LS_interp and LS_multi

% Evaluate the Objective and Gradient at the Initial Step
if nargout == 6
    [f_new,g_new,H] = funObj(x + t*d,varargin{:});
else
    [f_new,g_new] = funObj(x+t*d,varargin{:});
end
funEvals = 1;

while f_new > fr + c1*t*gtd || ~isLegal(f_new)
    temp = t;

    if LS_interp == 0 || ~isLegal(f_new)
        % Ignore value of new point
        if debug
            fprintf('Fixed BT\n');
        end
        t = 0.5*t;
    elseif LS_interp == 1 || ~isLegal(g_new)
        % Use function value at new point, but not its derivative
        if funEvals < 2 || LS_multi == 0 || ~isLegal(f_prev)
            % Backtracking w/ quadratic interpolation based on two points
            if debug
                fprintf('Quad BT\n');
            end
            t = polyinterp([0 f gtd; t f_new sqrt(-1)],doPlot,0,t);
        else
            % Backtracking w/ cubic interpolation based on three points
            if debug
                fprintf('Cubic BT\n');
            end
            t = polyinterp([0 f gtd; t f_new sqrt(-1); t_prev f_prev sqrt(-1)],doPlot,0,t);
        end
    else
        % Use function value and derivative at new point

        if funEvals < 2 || LS_multi == 0 || ~isLegal(f_prev)
            % Backtracking w/ cubic interpolation w/ derivative
            if debug
                fprintf('Grad-Cubic BT\n');
            end
            t = polyinterp([0 f gtd; t f_new g_new'*d],doPlot,0,t);
        elseif ~isLegal(g_prev)
            % Backtracking w/ quartic interpolation 3 points and derivative
            % of two
            if debug
                fprintf('Grad-Quartic BT\n');
            end
            t = polyinterp([0 f gtd; t f_new g_new'*d; t_prev f_prev sqrt(-1)],doPlot,0,t);
        else
            % Backtracking w/ quintic interpolation of 3 points and derivative
            % of two
            if debug
                fprintf('Grad-Quintic BT\n');
            end
            t = polyinterp([0 f gtd; t f_new g_new'*d; t_prev f_prev g_prev'*d],doPlot,0,t);
         end
    end

    % Adjust if change in t is too small/large
    if t < temp*1e-3
        if debug
            fprintf('Interpolated Value Too Small, Adjusting\n');
        end
        t = temp*1e-3;
    elseif t > temp*0.6
        if debug
            fprintf('Interpolated Value Too Large, Adjusting\n');
        end
        t = temp*0.6;
    end

    % Store old point if doing three-point interpolation
    if LS_multi
        f_prev = f_new;
        t_prev = temp;
        if LS_interp == 2
            g_prev = g_new;
        end
    end

    if ~saveHessianComp && nargout == 6
        [f_new,g_new,H] = funObj(x + t*d,varargin{:});
    else
        [f_new,g_new] = funObj(x + t*d,varargin{:});
    end
    funEvals = funEvals+1;

    % Check whether step size has become too small
    if max(abs(t*d)) <= progTol
        if debug
            fprintf('Backtracking Line Search Failed\n');
        end
        t = 0;
        f_new = f;
        g_new = g;
        break;
    end
end

% Evaluate Hessian at new point
if nargout == 6 && funEvals > 1 && saveHessianComp
    [f_new,g_new,H] = funObj(x + t*d,varargin{:});
    funEvals = funEvals+1;
end

x_new = x + t*d;

end

author: zouchuhang

repo: LayoutNet


text:
function theta = polygon3dNormalAngle(points, ind)
%POLYGON3DNORMALANGLE Normal angle at a vertex of the 3D polygon.
%
%   THETA = polygon3DNormalAngle(POLYGON, IND)
%   where POLYGON is a set of points, and IND is index of a point in
%   polygon. The function compute the angle of the normal cone localized at
%   this vertex.
%   If IND is a vector of indices, normal angle is computed for each vertex
%   specified by IND.
%
%   Example
%   % create an equilateral triangle in space
%   poly3d = [1 1 0;-1 0 1;0 -1 -1];
%   % compute each normal angle
%   theta = polygon3dNormalAngle(poly3d, 1:size(poly3d, 1));
%   % sum of normal angles must be equal to 2*PI for simple polygons
%   sum(theta)
%
%   IMPORTANT NOTE: works only for convex angles ! ! ! !
%
%   See also
%   polygons3d, faceNormalAngle
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2005-11-30
% Copyright 2005 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

% number of points
np = size(points, 1);

% number of angles to compute
nv = length(ind);

theta = zeros(nv, 1);

for i=1:nv
    p0 = points(ind(i), :);

    if ind(i)==1
        p1 = points(np, :);
    else
        p1 = points(ind(i)-1, :);
    end

    if ind(i)==np
        p2 = points(1, :);
    else
        p2 = points(ind(i)+1, :);
    end

    theta(i) = pi - anglePoints3d(p1, p0, p2);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function binary( varargin )

%BINARY Declares an binary variable.
%   BINARY VARIABLE x
%   where x is a valid MATLAB variable name, declares a scalar variable for
%   the current CVX problem, and constrains it to the set {0,1}. Put
%   another way, it is equivalent to
%       integer variable x
%       0 <= x <= 1
%
%   BINARY VARIABLE x(SZ), where SZ is a size vector, declares an array of
%   size SZ and constrains each element to be binary. Structure modifiers
%   such as "symmetric", "toeplitz", etc. are also permitted.
%
%   BINARY VARIABLES x y(SZ) z ... can be used to declare multiple binary
%   variables. Note however that structure modifiers are not permitted in
%   this case.
%
%   Obviously, binary variables are NOT convex. Problems with binary
%   variables can only be handled by solvers with explicit support for
%   them; in particular, neither of the free solvers SeDuMi nor SDPT3
%   provide binary support.
%
%   Examples:
%      binary variable x
%      binary variable x(100)
%
%   See also INTEGER, VARIABLE, VARIABLES.

if nargin < 2,
    error( 'Incorrect syntax for BINARY VARIABLE(S). Type HELP BINARY for details.' );
elseif ~iscellstr( varargin ),
    error( 'All arguments must be strings.' );
elseif strcmp( varargin{1}, 'variable' ),
    evalin( 'caller', sprintf( '%s ', 'variable', varargin{2:end}, 'binary' ) );
elseif strcmp( varargin{1}, 'variables' ),
    for k = 2 : nargin,
        evalin( 'caller', sprintf( '%s ', 'variable', varargin{k}, 'binary' ) );
    end
else
    error( 'Incorrect syntax for BINARY VARIABLE(S). Type HELP BINARY for details.' );
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [vol, labels] = imVolume(img, varargin)
% Volume of regions within a 3D binary or label image.
%
%   V = imVolume(IMG);
%   Compute volume of the image. IMG is either a binary image, or a label
%   image. In the case of a label image, the area of each labeled area is
%   returned in a column vector with as many elements as the number of
%   labels.
%
%   V = imVolume(IMG, SCALE);
%   Also specify scale of image tile. SCALE si a 1-by-3 array, containing
%   voxel size in each direction.
%
%   See Also
%   imVolumeDensity, imSurfaceArea, imMeanBreadth
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2010-01-15,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRAE - Cepia Software Platform.

%% Parse input arguments

% check image dimension
if ndims(img) ~= 3
    error('first argument should be a 3D binary or label image');
end

% the labels to compute
labels = [];

% default spatial calibration
delta = [1 1 1];

% parse parameter name-value pairs
while ~isempty(varargin)
    var1 = varargin{1};
    varargin(1) = [];

    if isnumeric(var1)
        % option can be number of directions, spatial calibration, or list
        % of labels
        if all(size(var1) == [1 3])
            % spatial calibration
            delta = var1;
        elseif ~isscalar(var1) && size(var1, 2) == 1
            % list of labels to compute
            labels = var1;
        else
            error('Unable to parse input argument');
        end
    else
        error('Expect numeric input only');
    end
end

%% Process label images

% in case of a label image, return a vector with a set of results
if ~islogical(img)
    % extract labels if necessary (considers 0 as background)
    if isempty(labels)
        labels = imFindLabels(img);
    end

    % allocate memory
    vol = zeros(length(labels), 1);

    % iterate over labels
    for i = 1:length(labels)
        vol(i) = imVolume(img==labels(i), delta);
    end
    return;
end

%% Process binary images

% in case of binary image, compute only one label...
labels = 1;

% compute area, multiplied by image resolution
vol = sum(img(:)) * prod(delta);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function pot = combine_pots(pot1, pot2)
% COMBINE_POTS combine two potentials
% pot = combine_pots(pot1, pot2)

% Reduce both potentials before trying to combine them.
% Cf. "Stable Local computation with Conditional Gaussian Distributions", page 9
% Consider again two potentials with minimal tail

% Guarantee minimal tails. If pot1 or pot2 are minimal, they are not changed
pot1 = reduce_pot(pot1);
pot2 = reduce_pot(pot2);

%if the intersect set of these two potentials' head conts. combination is undifined
if ~isempty( myintersect(pot1.cheaddom, pot2.cheaddom) )
    return;
end

if  isempty( myintersect(pot1.domain, pot2.cheaddom) ) | isempty( myintersect(pot2.domain, pot1.cheaddom))
    % if satisfy the condition of directed combine
    pot = direct_combine_pots(pot1, pot2);
else
    % perform recursive combine
    pot = recursive_combine_pots(pot1, pot2);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function physio = tapas_physio_new_philips(physio)
% Initializes physio-properties for Philips 3T Achieva system with good
% ECG-data
%
%   physio = tapas_physio_new_philips(physio)
%
% IN
%
% OUT
%
% EXAMPLE
%   tapas_physio_new_philips
%
%   See also

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

log_files   = physio.log_files;
scan_timing = physio.scan_timing;
model       = physio.model;
preproc     = physio.preproc;
verbose     = physio.verbose;

sqpar       = scan_timing.sqpar;
sync        = scan_timing.sync;

model.type = 'RETROICOR';
model.order = struct('c',3,'r',4,'cr',1, 'orthogonalise', 'none');
log_files.sampling_interval = 1/496; % for WiFi 1/496, cable: 1/500

sqpar.Ndummies = 0; % PARAM
sqpar.Nscans = 160; % PARAM
sqpar.Nprep = [];
sqpar.time_slice_to_slice = [];

model.order.c = 3;
model.order.r = 4;
model.order.cr = 1;
model.type = 'RETROICOR';

verbose.level = 2; % PARAM
verbose.fig_output_file ='physio_output.png';
verbose.use_tabs = 0;

sync.method = 'gradient_log';
sync.grad_direction = 'z';
sync.slice = 1500;
sync.zero = 1400;
sync.vol_spacing = [];

preproc.cardiac.modality = 'ECG';
preproc.cardiac.initial_cpulse_select.method = 'load_from_logfile'; 'auto_matched';
preproc.cardiac.posthoc_cpulse_select.method = 'off';

scan_timing.sqpar = sqpar;
scan_timing.sync = sync;

physio.log_files = log_files;
physio.scan_timing = scan_timing;
physio.preproc = preproc;
physio.model = model;
physio.verbose = verbose;

author: translationalneuromodeling

repo: tapas


text:
function imdb = setupVoc(datasetDir, varargin)
% SETUPVOC    Setup PASCAL VOC data
%   IMDB = SETUPVOC(DATASETDIR, 'EDITION', '2007') setups the
%   PASCAL VOC 2007 data. This is similar to SETUPGENERIC(), but adapted
%   to the peculiarities of this dataset. In particular, the
%   difficult image flag and the fact that multiple labels apply to
%   each image are supported.
%
%   Note that only the PASCAL VOC 2007 data comes with test images
%   and labels. For the other editions, setting up the test images
%   cannot be automatized due to restrictions in the distribution.
%
%   See also: SETUPGENERIC().

% Author: Andrea Vedaldi

% Copyright (C) 2013 Andrea Vedaldi
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.edition = '2007' ;
opts.autoDownload = true ;
opts.lite = false ;
opts = vl_argparse(opts, varargin) ;

switch opts.edition
  case '2007'
    urls = {'http://pascallin.ecs.soton.ac.uk/challenges/VOC/voc2007/VOCtrainval_06-Nov-2007.tar', ...
            'http://pascallin.ecs.soton.ac.uk/challenges/VOC/voc2007/VOCtest_06-Nov-2007.tar'} ;
  case '2008'
    urls = {'http://pascallin.ecs.soton.ac.uk/challenges/VOC/voc2008/VOCtrainval_14-Jul-2008.tar'} ;
  case '2009'
    urls = {'http://pascallin.ecs.soton.ac.uk/challenges/VOC/voc2009/VOCtrainval_11-May-2009.tar'} ;
  case '2010'
    urls = {'http://pascallin.ecs.soton.ac.uk/challenges/VOC/voc2010/VOCtrainval_03-May-2010.tar'} ;
  case '2011'
    urls = {'http://pascallin.ecs.soton.ac.uk/challenges/VOC/voc2011/VOCtrainval_25-May-2011.tar'} ;
  case '2012'
    urls = {'http://pascallin.ecs.soton.ac.uk/challenges/VOC/voc2012/VOCtrainval_11-May-2012.tar'} ;
  otherwise
    error('Unknown PASCAL VOC edition ''%s''.', opts.edition) ;
end

% Download and unpack
vl_xmkdir(datasetDir) ;
if exist(fullfile(datasetDir, ['VOC' opts.edition], 'Annotations'))
  % ok
elseif exist(fullfile(datasetDir, 'VOCdevkit', ['VOC' opts.edition], 'Annotations'))
  % ok
  datasetDir = fullfile(datasetDir, 'VOCdevkit') ;
elseif opts.autoDownload
  for i = 1:length(urls)
    fprintf('Downloading VOC data ''%s'' to ''%s''. This will take a while.', urls{i}, datasetDir) ;
    untar(urls{i}, datasetDir) ;
  end
  datasetDir = fullfile(datasetDir, 'VOCdevkit') ;
else
  error('VOC data not found in %s', datasetDir) ;
end

imdb.images.id = [] ;
imdb.images.set = uint8([]) ;
imdb.images.name = {} ;
imdb.images.size = zeros(2,0) ;
imdb.meta.sets = {'train', 'val', 'test'} ;
imdb.meta.classes = {'aeroplane', 'bicycle', 'bird', 'boat', 'bottle', 'bus', 'car', ...
  'cat', 'chair', 'cow', 'diningtable', 'dog', 'horse', 'motorbike', ...
  'person', 'pottedplant', 'sheep', 'sofa', 'train', 'tvmonitor'} ;
imdb.imageDir = fullfile(datasetDir, ['VOC', opts.edition], 'JPEGImages') ;

% Get the list of images
map = containers.Map() ;
j = 0 ;
for si = 1:numel(imdb.meta.sets)
  setName = imdb.meta.sets{si} ;
  annoPath = fullfile(datasetDir, ['VOC' opts.edition], ...
    'ImageSets', 'Main', ...
    [setName '.txt']) ;
  fprintf('%s: reading %s\n', mfilename, annoPath) ;
  [names,labels] = textread(annoPath, '%s %f') ;
  for i=1:length(names)
    j = j + 1 ;
    map(names{i}) = j ;
    imdb.images.id(j) = j ;
    imdb.images.set(j) = si ;
    imdb.images.name{j} = [names{i} '.jpg'] ;
    %info = imfinfo(fullfile(imdb.imageDir, imdb.images.name{j})) ;
    %imdb.images.size(:,j) = [info.Width ; info.Height] ;
    %fprintf('\radded %s', imdb.images.name{j}) ;
  end
end

% Get class labels
for ci = 1:length(imdb.meta.classes)
  imdb.classes.imageIds{ci} = [] ;
  imdb.classes.difficult{ci} = false(0) ;
end
for si = 1:numel(imdb.meta.sets)
  for ci = 1:length(imdb.meta.classes)
    setName = imdb.meta.sets{si} ;
    className = imdb.meta.classes{ci} ;
    annoPath = fullfile(datasetDir, ['VOC' opts.edition], ...
      'ImageSets', 'Main', ...
      [className '_' setName '.txt']) ;
    fprintf('%s: reading %s\n', mfilename, annoPath) ;
    [names,labels] = textread(annoPath, '%s %f') ;
    for i = 1:numel(names)
      j = map(names{i}) ;
      if labels(i) >= 0
        imdb.classes.imageIds{ci}(end+1) = j ;
        imdb.classes.difficult{ci}(end+1) = (labels(i) == 0) ;
      end
    end
  end
end

if opts.lite
  ok = {} ;
  for c = 1:3
    trainIds = intersect(imdb.images.id(imdb.images.set == 1), imdb.classes.imageIds{c}) ;
    testIds = intersect(imdb.images.id(imdb.images.set == 3), imdb.classes.imageIds{c}) ;

    ok{end+1} = vl_colsubset(find(ismember(imdb.images.id, trainIds)), 5) ;
    ok{end+1} = vl_colsubset(find(ismember(imdb.images.id, testIds)), 5) ;
  end
  ok = unique(cat(2, ok{:})) ;
  imdb.meta.classes = imdb.meta.classes(1:3) ;
  imdb.classes.imageIds = imdb.classes.imageIds(1:3) ;
  imdb.classes.difficult = imdb.classes.difficult(1:3) ;
  imdb.images.id = imdb.images.id(ok) ;
  imdb.images.name = imdb.images.name(ok) ;
  imdb.images.set = imdb.images.set(ok) ;
  for c = 1:3
    ok = ismember(imdb.classes.imageIds{c}, imdb.images.id) ;
    imdb.classes.imageIds{c} = imdb.classes.imageIds{c}(ok) ;
    imdb.classes.difficult{c} = imdb.classes.difficult{c}(ok) ;
  end
end

author: jianxiongxiao

repo: ProfXkit


text:
%%*******************************************************************
%%  Read in a problem in SDPA sparse format.
%%
%%  [blk,At,C,b] = read_sdpa(fname)
%%
%%  Input: fname = name of the file containing SDP data in
%%                 SDPA foramt.
%%  Important: the data is assumed to contain only
%%             semidefinite and linear blocks.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [blk,At,C,b] = read_sdpa(fname,deblksize)

   if (nargin<2); deblksize = 100; end
%%
%%  Open the file for input
%%
   compressed = 0;
   if exist(fname)
      fid = fopen(fname,'r');
   elseif exist([fname,'.Z']);
      compressed = 1;
      unix(['uncompress ',fname,'.Z']);
      fid = fopen(fname,'r');
   elseif exist([fname,'.gz']);
      compressed = 2;
      unix(['gunzip ',fname,'.gz']);
      fid = fopen(fname,'r');
   else
      fprintf('*** Problem not found, please specify the correct path or problem name. \n');
      blk = []; At = []; C = []; b = [];
      return;
   end
%%
%%  Clean up special characters and comments from the file
%%
   [datavec,count] = fscanf(fid,'%c');
   linefeeds = findstr(datavec,char(10));
   comment_chars = '*"=';
   cumidx = [];
   for i=1:length(comment_chars)
      idx = findstr(datavec,comment_chars(i));
      cumidx = [cumidx,idx];
   end
   for j=length(cumidx):-1:1
      if (cumidx(j)==1) | (strcmp(datavec(cumidx(j)-1),char(10)))
         datavec(cumidx(j):linefeeds(min(find(cumidx(j)<linefeeds))))='';
      else
         datavec(cumidx(j):linefeeds(min(find(cumidx(j)<linefeeds)))-1)='';
      end
   end
   special_chars = ',{}()';
   cumidx=[];
   for i=1:length(special_chars)
      idx = findstr(datavec,special_chars(i));
      cumidx = [cumidx,idx];
   end
   datavec(cumidx) = blanks(length(cumidx));
   clear linefeeds;
%%
%% Close the file
%%
   fclose('all');
   if compressed==1; unix(['compress ',fname]); end;
   if compressed==2; unix(['gzip ',fname]); end;
%%
%%  Next, read in basic problem size parameters.
%%
   datavec = sscanf(datavec,'%f');
   if size(datavec,1) < size(datavec,2); datavec = datavec'; end;
   m = datavec(1);
   numblk  = datavec(2);
   blksize = datavec(2+[1:numblk]);
   if size(blksize,1) > size(blksize,2); blksize = blksize'; end
%%
%% Get input  b.
%%
   idxstrt = 2+numblk;
   b = datavec(idxstrt+[1:m]);
   idxstrt = idxstrt+m;
   b = -b;
%%
%% Construct blk
%%
   spblkidxtmp = find( (blksize>1) & (blksize < deblksize) );
   spblkidxtmp = sort(spblkidxtmp);
   deblkidx = find( (blksize<=1) | (blksize >= deblksize) );
   denumblk = length(deblkidx);
   linblkidx = zeros(1,denumblk);
   for p = 1:denumblk
      n = blksize(deblkidx(p));
      if (n > 1);
         blk{p,1} = 's'; blk{p,2} = n;
         n2 = n*(n+1)/2;
         At{p,1} = sparse(n2,m);
         C{p,1} = sparse(n,n);
      else
         linblkidx(p) = p;
         blk{p,1} = 'l'; blk{p,2} = abs(n);
         At{p,1} = sparse(abs(n),m);
         C{p,1} = sparse(abs(n),1);
      end
   end
   if ~isempty(spblkidxtmp)
      maxnumblk = 200;
      spnumblk = ceil(length(spblkidxtmp)/maxnumblk);
      for q = 1:spnumblk
         if (q < spnumblk)
            spblkidxall{q} = spblkidxtmp([(q-1)*maxnumblk+1: q*maxnumblk]);
         else
            spblkidxall{q} = spblkidxtmp([(q-1)*maxnumblk+1: length(spblkidxtmp)]);
         end
         tmp = blksize(spblkidxall{q});
         blk{denumblk+q,1} = 's';
         blk{denumblk+q,2} = tmp;
         n2 = sum(tmp.*(tmp+1))/2;
         At{denumblk+q,1} = sparse(n2,m);
         C{denumblk+q,1} = sparse(sum(tmp),sum(tmp));
      end
   else
      spnumblk = 0;
   end
   linblkidx(denumblk+[1:spnumblk]) = zeros(1,spnumblk);
%%
%% Construct single blocks of A,C
%%
   len = length(datavec);
   Y = reshape(datavec(idxstrt+1:len),5,(len-idxstrt)/5)';
   clear datavec;
   Y = sortrows(Y,[1 2]);
   matidx = [0; find(diff(Y(:,1)) ~= 0); size(Y,1)];
%%
   for k = 1:length(matidx)-1
      idx = [matidx(k)+1 : matidx(k+1)];
      Ytmp  = Y(idx,1:5);
      matno = Ytmp(1,1);
      Ytmp2 = Ytmp(:,2);
      for p = 1:denumblk
         n  = blksize(deblkidx(p));
         idx = find(Ytmp2 == deblkidx(p));
         ii = Ytmp(idx,3); jj = Ytmp(idx,4); vv =Ytmp(idx,5);
         len = length(idx);
         if (n > 1)
            idxtmp = find(ii > jj);
            if ~isempty(idxtmp);
               tmp = jj(idxtmp);
               jj(idxtmp) = ii(idxtmp); ii(idxtmp) = tmp;
            end
            tmp = -sparse(ii,jj,vv,n,n);
            tmp = tmp + triu(tmp,1)';
         else
            tmp = -sparse(ii,ones(len,1),vv,abs(n),1);
         end
         if (matno == 0)
            C{p,1} = tmp;
         else
            if (n > 1)
               At{p,1}(:,matno) = svec(blk(p,:),tmp,1);
            else
               At{p,1}(:,matno) = tmp;
            end
         end
      end
   end
%%
%% Construct big sparse block of A,C
%%
if (spnumblk > 0)
   Y1 = Y(:,1);
   diffY1 = find(diff([-1; Y1; inf]));
   for kk = 1:length(diffY1)-1
      idx = [diffY1(kk) : diffY1(kk+1)-1];
      matno = Y1(diffY1(kk));
      Ytmp  = Y(idx,1:5);
      Ytmp2 = Ytmp(:,2);
      maxYtmp2  = Ytmp2(length(Ytmp2));
      minYtmp2  = Ytmp2(1);
      diffYtmp2 = Ytmp2(find(diff([-1; Ytmp2])));
      for q = 1:spnumblk
         spblkidx    = spblkidxall{q};
         maxspblkidx = spblkidx(length(spblkidx));
         minspblkidx = spblkidx(1);
         count = 0;
         if (minYtmp2 <= maxspblkidx) & (maxYtmp2 >= minspblkidx)
            tmpblksize = blksize(spblkidx);
            n = sum(tmpblksize);
            cumspblksize = [0 cumsum(tmpblksize)];
            n2 = sum(tmpblksize.*(tmpblksize+1))/2;
            idx = zeros(n2,1); offset = zeros(n2,1);
            for t = [1:length(diffYtmp2)]
  	       p = find(spblkidx == diffYtmp2(t));
               if ~isempty(p)
	          idxtmp = find(Ytmp2 == spblkidx(p));
                  len = length(idxtmp);
     	          idx(count+[1:len]) = idxtmp;
                  offset(count+[1:len]) = cumspblksize(p);
                  count = count + len;
               end
            end
            idx = idx(1:count); offset = offset(1:count);
            ii = Ytmp(idx,3)+offset; jj = Ytmp(idx,4)+offset; vv = Ytmp(idx,5);
            idxtmp = find(ii > jj);
            if ~isempty(idxtmp);
               tmp = jj(idxtmp);
               jj(idxtmp) = ii(idxtmp); ii(idxtmp) = tmp;
            end
	    idxeq = find(ii==jj);
            tmp = spconvert([ii jj -vv; jj ii -vv; n n 0]) ...
                  + spconvert([ii(idxeq) jj(idxeq) vv(idxeq); n n 0]);
            if (matno == 0)
               C{denumblk+q,1} = tmp;
            else
               At{denumblk+q,1}(:,matno) = svec(blk(denumblk+q,:),tmp,1);
            end
         end
      end
   end
end
%%
%% put all linear blocks together as a single linear block
%%
   idx = find(linblkidx);
   if (length(idx) > 1)
      sdpidx = find(linblkidx==0);
      blktmp = 0; Atmp = []; Ctmp = [];
      for k = 1:length(idx)
          tmp = linblkidx(idx(k));
          blktmp = blktmp+blk{tmp,2};
          Atmp = [Atmp; At{tmp}];
          Ctmp = [Ctmp; C{tmp}];
      end
      At = At(sdpidx); C = C(sdpidx); blk = blk(sdpidx,:);
      len = length(sdpidx);
      blk(2:len+1,:) = blk;
      blk{1,1} = 'l'; blk{1,2} = blktmp;
      At(2:len+1,1) = At; C(2:len+1,1) = C;
      At{1,1} = Atmp; C{1,1} = Ctmp;
   end
%%******************************************************************

author: shamilmamedov


text:
% pop_newtimef() - Returns estimates and plots of event-related (log) spectral
%           perturbation (ERSP) and inter-trial coherence (ITC) phenomena
%           timelocked to a set of single-channel input epochs
%
% Usage:
%   >> pop_newtimef(EEG, typeplot);          % pop_up window
%   >> pop_newtimef(EEG, typeplot, lastcom); % pop_up window
%   >> pop_newtimef(EEG, typeplot, channel); % do not pop-up window
%   >> pop_newtimef(EEG, typeproc, num, tlimits,cycles,
%                        'key1',value1,'key2',value2, ... );
%
% Graphical interface:
%   "Channel/component number" - [edit box] this is the index of the data
%              channel or the index of the component for which to plot the
%              time-frequency decomposition.
%   "Sub-epoch time limits" - [edit box] sub epochs may be extracted (note that
%              this function aims at plotting data epochs not continuous data).
%              You may select the new epoch limits in this edit box.
%   "Use n time points" - [muliple choice list] this is the number of time
%              points to use for the time-frequency decomposition. The more
%              time points, the longer the time-frequency decomposition
%              takes to compute.
%   "Frequency limits" - [edit box] these are the lower and upper
%              frequency limit of the time-frequency decomposition. Instead
%              of limits, you may also enter a sequence of frequencies. For
%              example to compute the time-frequency decomposition at all
%              frequency between 5 and 50 hertz with 1 Hz increment, enter "1:50"
%   "Use limits, padding n" - [muliple choice list] "using limits" means
%              to use the upper and lower limits in "Frequency limits" with
%              a specific padding ratio (padratio argument of newtimef).
%              The last option "use actual frequencies" forces newtimef to
%              ignore the padratio argument and use the vector of frequencies
%              given as input in the "Frequency limits" edit box.
%   "Log spaced" - [checkbox] you may check this box to compute log-spaced
%              frequencies. Note that this is only relevant if you specify
%              frequency limits (in case you specify actual frequencies,
%              this parameter is ignored).
%   "Use divisive baseline" - [muliple choice list] there are two types of
%              baseline correction, additive (the baseline is subtracted)
%              or divisive (the data is divided by the baseline values).
%              The choice is yours. There is also the option to perform
%              baseline correction in single trials. See the 'trialbase' "full"
%              option in the newtimef.m documentation for more information.
%   "No baseline" - [checkbox] check this box to compute the raw time-frequency
%              decomposition with no baseline removal.
%   "Wavelet cycles" - [edit box] specify the number of cycle at the lowest
%              and highest frequency. Instead of specifying the number of cycle
%              at the highest frequency, you may also specify a wavelet
%              "factor" (see newtimef help message). In addition, it is
%              possible to specify actual wavelet cycles for each frequency
%              by entering a sequence of numbers.
%   "Use FFT" - [checkbox] check this checkbox to use FFT instead of
%              wavelet decomposition.
%   "ERSP color limits" - [edit box] set the upper and lower limit for the
%              ERSP image.
%   "see log power" - [checkbox] the log power values (in dB) are plotted.
%              Uncheck this box to plot the absolute power values.
%   "ITC color limits" - [edit box] set the upper and lower limit for the
%              ITC image.
%   "plot ITC phase" - [checkbox] check this box plot plot (overlayed on
%              the ITC amplitude) the polarity of the ITC complex value.
%   "Bootstrap significance level" - [edit box] use this edit box to enter
%              the p-value threshold for masking both the ERSP and the ITC
%              image for significance (masked values appear as light green)
%   "FDR correct" - [checkbox] this correct the p-value for multiple comparisons
%              (accross all time and frequencies) using the False Discovery
%              Rate method. See the fdr.m function for more details.
%   "Optional newtimef arguments" - [edit box] addition argument for the
%              newtimef function may be entered here in the 'key', value
%              format.
%   "Plot Event Related Spectral Power" - [checkbox] plot the ERSP image
%              showing event related spectral stimulus induced changes
%   "Plot Inter Trial Coherence" - [checkbox] plot the ITC image.
%   "Plot Curve at each frequency" - [checkbox] instead of plotting images,
%              it is also possible to display curves at each frequency.
%              This functionality is beta and might not work in all cases.
%
% Inputs:
%   INEEG    - input EEG dataset
%   typeproc - type of processing: 1 process the raw channel data
%                                  0 process the ICA component data
%   num      - component or channel number
%   tlimits  - [mintime maxtime] (ms) sub-epoch time limits to plot
%   cycles   -  > 0 --> Number of cycles in each analysis wavelet
%               = 0 --> Use FFTs (with constant window length
%                       at all frequencies)
%
% Optional inputs:
%    See the newtimef() function.
%
% Outputs: Same as newtimef(); no outputs are returned when a
%          window pops-up to ask for additional arguments
%
% Saving the ERSP and ITC output values:
%    Simply look up the history using the eegh function (type eegh).
%    Then copy and paste the pop_newtimef() command and add output args.
%    See the newtimef() function for a list of outputs. For instance,
% >> [ersp itc powbase times frequencies] = pop_newtimef( EEG, ....);
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: newtimef(), eeglab()

% Copyright (C) 2002 University of California San Diego
%
% 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% 01-25-02 reformated help & license -ad
% 03-08-02 add eeglab option & optimize variable sizes -ad
% 03-10-02 change newtimef call -ad
% 03-18-02 added title -ad & sm
% 04-04-02 added outputs -ad & sm

function varargout = pop_newtimef(EEG, typeproc, num, tlimits, cycles, varargin );

varargout{1} = '';
% display help if not enough arguments
% ------------------------------------
if nargin < 2
	help pop_newtimef;
	return;
end;
lastcom = [];
if nargin < 3
	popup = 1;
else
	popup = isstr(num) | isempty(num);
	if isstr(num)
		lastcom = num;
	end;
end;

% pop up window
% -------------
if popup
	[txt vars] = gethelpvar('newtimef.m');

    g = [1 0.3 0.6 0.4];
	geometry = { g g g g g g g g [0.975 1.27] [1] [1.2 1 1.2]};
    uilist = { ...
               { 'Style', 'text', 'string', fastif(typeproc, 'Channel number', 'Component number'), 'fontweight', 'bold'  } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,3,[],'1') 'tag' 'chan'} {} {} ...
               ...
			   { 'Style', 'text', 'string', 'Sub epoch time limits [min max] (msec)', 'fontweight', 'bold' } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,4,[],[ int2str(EEG.xmin*1000) ' ' int2str(EEG.xmax*1000) ]) 'tag' 'tlimits' } ...
               { 'Style', 'popupmenu', 'string', 'Use 50 time points|Use 100 time points|Use 150 time points|Use 200 time points|Use 300 time points|Use 400 time points' 'tag' 'ntimesout' 'value' 4} { } ...
			   ...
			   { 'Style', 'text', 'string', 'Frequency limits [min max] (Hz) or sequence', 'fontweight', 'bold' } ...
			   { 'Style', 'edit', 'string', '' 'tag' 'freqs'  } ...
               { 'Style', 'popupmenu', 'string', 'Use limits, padding 1|Use limits, padding 2|Use limits, padding 4|Use actual freqs.' 'tag' 'nfreqs' }  ...
               { 'Style', 'checkbox', 'string' 'Log spaced' 'value' 0 'tag' 'freqscale' } ...
			   ...
			   { 'Style', 'text', 'string', 'Baseline limits [min max] (msec) (0->pre-stim.)', 'fontweight', 'bold' } ...
			   { 'Style', 'edit', 'string', '0' 'tag' 'baseline' } ...
			   { 'Style', 'popupmenu',  'string', 'Use divisive baseline (DIV)|Use standard deviation (STD)|Use single trial DIV baseline|Use single trial STD baseline' 'tag' 'basenorm' } ...
               { 'Style', 'checkbox', 'string' 'No baseline' 'tag' 'nobase' } ...
               ...
               { 'Style', 'text', 'string', 'Wavelet cycles [min max/fact] or sequence', 'fontweight', 'bold' } ...
               { 'Style', 'edit', 'string', getkeyval(lastcom,5,[],'3 0.5') 'tag' 'cycle' } ...
               { 'Style', 'checkbox', 'string' 'Use FFT' 'value' 0 'tag' 'fft' } ...
               { } ...
			   ...
			   { 'Style', 'text', 'string', 'ERSP color limits [max] (min=-max)', 'fontweight', 'bold' } ...
               { 'Style', 'edit', 'string', '' 'tag' 'erspmax'} ...
               { 'Style', 'checkbox', 'string' 'see log power (set)' 'tag' 'scale' 'value' 1} {} ...
			   ...
			   { 'Style', 'text', 'string', 'ITC color limits [max]', 'fontweight', 'bold' } ...
               { 'Style', 'edit', 'string', '' 'tag' 'itcmax'} ...
               { 'Style', 'checkbox', 'string' 'plot ITC phase (set)' 'tag' 'plotphase' } {} ...
			   ...
			   { 'Style', 'text', 'string', 'Bootstrap significance level (Ex: 0.01 -> 1%)', 'fontweight', 'bold' } ...
               { 'Style', 'edit', 'string', getkeyval(lastcom,'alpha') 'tag' 'alpha'} ...
               { 'Style', 'checkbox', 'string' 'FDR correct (set)' 'tag' 'fdr' } {} ...
			   ...
			   { 'Style', 'text', 'string', 'Optional newtimef() arguments (see Help)', 'fontweight', 'bold', ...
				 'tooltipstring', 'See newtimef() help via the Help button on the right...' } ...
			   { 'Style', 'edit', 'string', '' 'tag' 'options' } ...
			   {} ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotersp','present',0), 'string', ...
				 'Plot Event Related Spectral Power', 'tooltipstring', ...
				 'Plot log spectral perturbation image in the upper panel' 'tag' 'plotersp' } ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotitc','present',0), 'string', ...
				 'Plot Inter Trial Coherence', 'tooltipstring', ...
				 'Plot the inter-trial coherence image in the lower panel' 'tag' 'plotitc' } ...
			   { 'Style', 'checkbox', 'value', 0, 'string', ...
				 'Plot curve at each frequency' 'tag' 'plotcurve' } ...
			 };
      % { 'Style', 'edit', 'string', '''padratio'', 4, ''plotphase'', ''off''' } ...
			   %{ 'Style', 'text', 'string',  '[set] -> Plot ITC phase sign', 'fontweight', 'bold', ...
			%	 'tooltipstring', ['Plot the sign (+/-) of inter-trial coherence phase' 10 ...
			%		'as red (+) or blue (-)'] } ...
			%   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotphase','present',1) } { } ...

	[ tmp1 tmp2 strhalt result ] = inputgui( geometry, uilist, 'pophelp(''pop_newtimef'');', ...
					   fastif(typeproc, 'Plot channel time frequency -- pop_newtimef()', ...
							  'Plot component time frequency -- pop_newtimef()'));
	if length( tmp1 ) == 0 return; end;

	if result.fft,      result.cycle = '0'; end;
	if result.nobase,   result.baseline = 'NaN'; end;

	num	     = eval( [ '[' result.chan    ']' ] );
	tlimits	 = eval( [ '[' result.tlimits ']' ] );
	cycles	 = eval( [ '[' result.cycle   ']' ] );
    freqs    = eval( [ '[' result.freqs   ']' ] );
    %result.ncycles == 2 is ignored

    % add topoplot
    % ------------
    options = [];
    if isfield(EEG.chanlocs, 'theta') && ~isempty(EEG.chanlocs(num).theta)
        if ~isfield(EEG, 'chaninfo'), EEG.chaninfo = []; end;
        if typeproc == 1
            if isempty(EEG.chanlocs), caption = [ 'Channel ' int2str(num) ]; else caption = EEG.chanlocs(num).labels; end;
            options = [options ', ''topovec'', ' int2str(num) ...
                        ', ''elocs'', EEG.chanlocs, ''chaninfo'', EEG.chaninfo, ''caption'', ''' caption '''' ];
        else
            options = [options ', ''topovec'', EEG.icawinv(:,' int2str(num) ...
                       '), ''elocs'', EEG.chanlocs, ''chaninfo'', EEG.chaninfo, ''caption'', [''IC ' num2str(num) ''']' ];
      end;
    end;

	if ~isempty( result.baseline ),  options = [ options ', ''baseline'',[' result.baseline ']' ]; end;
    if ~isempty( result.alpha ),     options = [ options ', ''alpha'',' result.alpha ];   end;
	if ~isempty( result.options ),   options = [ options ',' result.options ];            end;
	if ~isempty( result.freqs ),     options = [ options ', ''freqs'', [' result.freqs ']'   ]; end;
	if ~isempty( result.erspmax ),   options = [ options ', ''erspmax'', [' result.erspmax ']' ]; end;
	if ~isempty( result.itcmax ),    options = [ options ', ''itcmax'','  result.itcmax ];      end;
	if ~result.plotersp,             options = [ options ', ''plotersp'', ''off''' ];     end;
	if ~result.plotitc,              options = [ options ', ''plotitc'' , ''off''' ];     end;
	if result.plotcurve,             options = [ options ', ''plottype'', ''curve''' ];   end;
	if result.fdr,                   options = [ options ', ''mcorrect'', ''fdr''' ];     end;
	if result.freqscale,             options = [ options ', ''freqscale'', ''log''' ];    end;
	if ~result.plotphase,            options = [ options ', ''plotphase'', ''off''' ];    end;
	if ~result.scale,                options = [ options ', ''scale'', ''abs''' ];        end;
    if result.ntimesout == 1,        options = [ options ', ''ntimesout'', 50' ];         end;
    if result.ntimesout == 2,        options = [ options ', ''ntimesout'', 100' ];        end;
    if result.ntimesout == 3,        options = [ options ', ''ntimesout'', 150' ];        end;
    if result.ntimesout == 5,        options = [ options ', ''ntimesout'', 300' ];        end;
    if result.ntimesout == 6,        options = [ options ', ''ntimesout'', 400' ];        end;
    if result.nfreqs == 1,           options = [ options ', ''padratio'', 1' ];           end;
    if result.nfreqs == 2,           options = [ options ', ''padratio'', 2' ];           end;
    if result.nfreqs == 3,           options = [ options ', ''padratio'', 4' ];           end;
    if result.nfreqs == 4,           options = [ options ', ''nfreqs'', ' int2str(length(freqs)) ]; end;
    if result.basenorm == 2,         options = [ options ', ''basenorm'', ''on''' ];      end;
    if result.basenorm == 4,         options = [ options ', ''basenorm'', ''on''' ];      end;
    if result.basenorm >= 3,         options = [ options ', ''trialbase'', ''full''' ];      end;

    % add title
    % ---------
	if isempty( findstr(  '''title''', result.options))
        if ~isempty(EEG.chanlocs) & typeproc
            chanlabel = EEG.chanlocs(num).labels;
        else
            chanlabel = int2str(num);
        end;
	end;

    % compute default winsize
    % -----------------------
    if EEG.xmin < 0 && isempty(findstr( '''winsize''', result.options)) && isempty( result.freqs )
        fprintf('Computing window size in pop_newtimef based on half of the length of the baseline period');
        options = [ options ', ''winsize'', ' int2str(-EEG.xmin*EEG.srate) ];
    end;

	figure; try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end;
else
    options = [ ',' vararg2str(varargin) ];
end;

% compute epoch limits
% --------------------
if isempty(tlimits)
	tlimits = [EEG.xmin, EEG.xmax]*1000;
end;
pointrange1 = round(max((tlimits(1)/1000-EEG.xmin)*EEG.srate, 1));
pointrange2 = round(min((tlimits(2)/1000-EEG.xmin)*EEG.srate, EEG.pnts));
pointrange = [pointrange1:pointrange2];

% call function sample either on raw data or ICA data
% ---------------------------------------------------
if typeproc == 1
	tmpsig = EEG.data(num,pointrange,:);
else
	if ~isempty( EEG.icasphere )
        if ~isempty(EEG.icaact)
    		tmpsig = EEG.icaact(num,pointrange,:);
 	    else
            tmpsig = (EEG.icaweights(num,:)*EEG.icasphere)*reshape(EEG.data(:,pointrange,:), EEG.nbchan, EEG.trials*length(pointrange));
        end;
	else
		error('You must run ICA first');
	end;
end;
tmpsig = reshape( tmpsig, length(num), size(tmpsig,2)*size(tmpsig,3));

% outputs
% -------
outstr = '';
if ~popup
    for io = 1:nargout, outstr = [outstr 'varargout{' int2str(io) '},' ]; end;
    if ~isempty(outstr), outstr = [ '[' outstr(1:end-1) '] =' ]; end;
end;

% plot the datas and generate output command
% --------------------------------------------
if length( options ) < 2
    options = '';
end;
if nargin < 4
    varargout{1} = sprintf('figure; pop_newtimef( %s, %d, %d, [%s], [%s] %s);', inputname(1), typeproc, num, ...
			int2str(tlimits), num2str(cycles), options);
end;
com = sprintf('%s newtimef( tmpsig(:, :), length(pointrange), [tlimits(1) tlimits(2)], EEG.srate, cycles %s);', outstr, options);
eval(com)

return;

% get contextual help
% -------------------
function txt = context(var, allvars, alltext);
	loc = strmatch( var, allvars);
	if ~isempty(loc)
		txt= alltext{loc(1)};
	else
		disp([ 'warning: variable ''' var ''' not found']);
		txt = '';
	end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function pass = test_parSimp(varargin)

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

S = loadtests;

pass = zeros(numel(S), 1);
for k = 1:numel(S)
    str = stringParser.parSimp(S{k}{1});
    pass(k) = strcmp(str,S{k}{2});
end

function S = loadtests
% BVPs
S{1} = {'((0.02.*diff(u,2)+diff(u))+u)',
    '0.02.*diff(u,2)+diff(u)+u'};
S{2} = {'((0.01.*diff(u,2)-x.*u)-1)',
    '0.01.*diff(u,2)-x.*u-1'};
S{3} = {'((x.^(2).*diff(u,2)+x.*diff(u))+(x.^(2)-3.^(2)).*u)',
    'x.^2.*diff(u,2)+x.*diff(u)+(x.^2-3.^2).*u'};
S{4} = {'(((0.01.*diff(u,2)+2.*(1-x.^(2)).*u)+u.^(2))-1)',
    '0.01.*diff(u,2)+2.*(1-x.^2).*u+u.^2-1'};
S{5} = {'((diff(u,2)+(1.2+sign((10-abs(x)))).*u)-1)',
    'diff(u,2)+(1.2+sign((10-abs(x)))).*u-1'};
S{6} = {'((diff(u,2)+u)-u.^(2))',
    'diff(u,2)+u-u.^2'};
S{7} = {'(diff(u,4)-(diff(u).*diff(u,2)-u.*diff(u,3)))',
    'diff(u,4)-diff(u).*diff(u,2)+u.*diff(u,3)'};
S{8} = {'(diff(u,2)+.87.*exp(u))',
    'diff(u,2)+.87.*exp(u)'};
S{9} = {'((.0005.*diff(u,2)+x.*(x.^(2)-0.5).*diff(u))+3.*(x.^(2)-0.5).*u)',
    '.0005.*diff(u,2)+x.*(x.^2-0.5).*diff(u)+3.*(x.^2-0.5).*u'};
S{10} = {'((.01.*diff(u,2)+u.*diff(u))-u)',
    '.01.*diff(u,2)+u.*diff(u)-u'};
S{11} = {'(diff(u,2)+sin(u))',
    'diff(u,2)+sin(u)'};
S{12} = {'((0.05.*diff(u,2)+diff(u).^(2))-1)',
    '0.05.*diff(u,2)+diff(u).^2-1'};
S{13} = {'(diff(u,2)-8.*sinh(8.*u))',
        'diff(u,2)-8.*sinh(8.*u)'};
S{14} = {'(diff(u,2)-(u-1).*(1+diff(u).^(2)).^(1.5))',
    'diff(u,2)-(u-1).*(1+diff(u).^2).^1.5'};
S{15} = {'((diff(u,2)-x.*sin(u))-1)',
    'diff(u,2)-x.*sin(u)-1'};
S{16} = {'((diff(u,2)-(1-u.^(2)).*diff(u))+u)',
    'diff(u,2)-(1-u.^2).*diff(u)+u'};
S{17} = {'(diff(u,2)-sin(v))',
    'diff(u,2)-sin(v)'};
S{18} = {'(diff(u,2)-sin(v))',
    'diff(u,2)-sin(v)'};
S{19} = {'(cos(u)+diff(v,2))',
    'cos(u)+diff(v,2)'};

% PDEs
S{20} = {'+(0.1.*diff(u,2)+diff(u))',
    '0.1.*diff(u,2)+diff(u)'};
S{21} = {'+((.01.*diff(u,2)+u)-u.^(3))',
    '.01.*diff(u,2)+u-u.^3'};
S{22} = {'+(-diff(u.^(2))+.02.*diff(u,2))',
    '-diff(u.^2)+.02.*diff(u,2)'};
S{23} = {'+(-.003.*diff(u,4)+diff((u.^(3)-u),2))',
    '-.003.*diff(u,4)+diff(u.^3-u,2)'};
S{24} = {'+0.1.*diff(u,2)',
    '0.1.*diff(u,2)'};
S{25} = {'+(.02.*diff(u,2)+cumsum(u).*sum(u))',
    '.02.*diff(u,2)+cumsum(u).*sum(u)'};
S{26} = {'+((u.*diff(u)-diff(u,2))-0.006.*diff(u,4))',
    'u.*diff(u)-diff(u,2)-0.006.*diff(u,4)'};

% Systems
S{27} = {'+((-u+(x+1).*v)+0.1.*diff(u,2))',
    '-u+(x+1).*v+0.1.*diff(u,2)'};
S{28} = {'+((u-(x+1).*v)+0.2.*diff(v,2))',
    'u-(x+1).*v+0.2.*diff(v,2)'};
S{29} = {'+(0.1.*diff(u,2)-100.*u.*v)',
    '0.1.*diff(u,2)-100.*u.*v'};
S{30} = {'+(.2.*diff(v,2)-100.*u.*v)',
    '.2.*diff(v,2)-100.*u.*v'};
S{31} = {'+(0.001.*diff(w,2)+200.*u.*v)',
    '0.001.*diff(w,2)+200.*u.*v'};
S{32} = {'+(diff(u,2)-v)',
    'diff(u,2)-v'};
S{33} = {'(diff(v,2)-u)',
    'diff(v,2)-u'};
S{34} = {'+(diff(u,2)+exp(-100.*x.^(2)).*sin(pi.*t))',
    'diff(u,2)+exp(-100.*x.^2).*sin(pi.*t)'};

% EIGs
S{35} = {'(diff(u,2)+diff(u))',
    'diff(u,2)+diff(u)'};
S{36} = {'(-diff(u,2)+1i.*x.^(2).*u)',
    '-diff(u,2)+1i.*x.^2.*u'};
S{37} = {'(-.1.*diff(u,2)+4.*(sign((x+1))-sign((x-1))).*u)',
    '-.1.*diff(u,2)+4.*(sign((x+1))-sign((x-1))).*u'};
S{38} = {'(-diff(u,2)+x.^(2).*u)',
    '-diff(u,2)+x.^2.*u'};
S{39} = {'(-diff(u,2)+5.*cos(2.*x).*u)',
    '-diff(u,2)+5.*cos(2.*x).*u'};
S{40} = {'((diff(u,2)+u.*x)+v)',
    'diff(u,2)+u.*x+v'};
S{41} = {'(diff(v,2)+sin(x).*u)',
    'diff(v,2)+sin(x).*u'};

% MISC

% Check exponential notation
S{42} = {'(diff(v,2)-1e-2*u)',
    'diff(v,2)-1e-2*u'};
S{43} = {'(diff(v,2)-1e+2*u)',
    'diff(v,2)-1e+2*u'};
S{44} = {'(diff(v,2)+1e-2*u)',
    'diff(v,2)+1e-2*u'};
S{45} = {'(diff(v,2)+1e+2*u)',
    'diff(v,2)+1e+2*u'};

% Check division (see GitHub #2357).
S{46} = {'diff(y)-(2/3)/4',
    'diff(y)-2/3/4'};
S{47} = {'diff(y)-2/(3/4)',
    'diff(y)-2/(3/4)'};
S{48} = {'diff(y)-1/(2*(y-1))',
    'diff(y)-1/(2*(y-1))'};

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function varargout = spm_diffeo(varargin)
% Mex function called for image registration stuff
%
%_______________________________________________________________________
%
% FORMAT u = spm_diffeo('vel2mom', v, param)
% v     - velocity (flow) field n1*n2*n3*3.
% param - 8 parameters (settings)
%         - [1][2][3] Voxel sizes
%         - [4][5][6][7][8] Regularisation parameters
%           - [4] Absolute displacements need to be penalised by a tiny
%                 amount.  The first element encodes the amount of
%                 penalty on these.  Ideally, absolute displacements
%                 should not be penalised, but it is usually necessary
%                 for technical reasons.
%           - [5] The `membrane energy' of the deformation is penalised,
%                 usually by a relatively small amount. This penalises
%                 the sum of squares of the derivatives of the velocity
%                 field (ie the sum of squares of the elements of the
%                 Jacobian tensors).
%           - [6] The `bending energy' is penalised (3rd element). This
%                 penalises the sum of squares of the 2nd derivatives of
%                 the velocity.
%           - [7][8] Linear elasticity regularisation is also included.
%                    The first parameter (mu) is similar to that for
%                    linear elasticity, except it penalises the sum of
%                    squares of the Jacobian tensors after they have been
%                    made symmetric (by averaging with the transpose).
%                    This term essentially penalises length changes,
%                    without penalising rotations.
%                    The final term also relates to linear elasticity,
%                    and is the weight that denotes how much to penalise
%                    changes to the divergence of the velocities (lambda).
%                    This divergence is a measure of the rate of volumetric
%                    expansion or contraction.
% u       - `momentum' field n1*n2*n3*3.
%
% Convert a velocity field to a momentum field by u = A*v, where
% A is the large sparse matrix encoding some form of regularisation.
% v and m are single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT v = spm_diffeo('mom2vel',g, param)
% v     - the solution n1*n2*n3*3
% g     - parameterisation of first derivatives
% param - 10 parameters (settings)
%         - [1][2][3] Voxel sizes
%         - [4][5][6][7][8] Regularisation settings (see vel2mom).
%         - [9] Number of Full Multigrid cycles.
%         - [10] Number of relaxation iterations per cycle.
%
% Solve equations using a Full Multigrid method.  See Press et al
% for more information.
% v = inv(A)*g
% g and v are both single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT v = spm_diffeo('fmg',H, g, param)
% v     - the solution n1*n2*n3*3
% H     - parameterisation of 2nd derivatives
% g     - parameterisation of first derivatives
% param - 10 parameters (settings)
%         - [1][2][3] Voxel sizes
%         - [4][5][6][7][8] Regularisation settings (see vel2mom).
%         - [9] Number of Full Multigrid cycles.
%         - [10] Number of relaxation iterations per cycle.
%
% Solve equations using a Full Multigrid method, but using Hessian of
% the matching term.  See Press et al for more information.
% v = inv(A+H)*g
% H, g and v are all single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT v = spm_diffeo('cgs',H, g, param)
% v     - the solution
% H     - parameterisation of 2nd derivatives
% g     - parameterisation of first derivatives
% param - 10 parameters (settings)
%         - [1][2][3] Voxel sizes
%         - [4][5][6][7][8] Regularisation settings (see vel2mom).
%         - [9] Tolerance.  Indicates required degree of accuracy.
%         - [10] Maximum number of iterations.
%
% This is for solving a set of equations using a conjugate gradient
% solver. This method is less efficient than the Full Multigrid, and
% is included for illustrative purposes.
% v = inv(A+H)*g
% H, g and v are all single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT F = spm_diffeo('kernel',d,prm)
% d   - image dimensions
% prm - 8 parameters (settings).
%       These are described above (for 'vel2mom').
% F   - The differential operator encoded as an image (or images).
%       Convolving a velocity field by this will give the momentum.
%
%_______________________________________________________________________
%
% FORMAT y3 = spm_diffeo('comp',y1,y2)
% y1, y2 - deformation fields n1*n2*n3*3.
% y3     - deformation field field n1*n2*n3*3.
%
% Composition of two deformations y3 = y1(y2)
% y1, y2 and y3 are single precision floating point.
%
%
% FORMAT [y3,J3] = spm_diffeo('comp', y1, y2, J1, J2)
% y1, y2 - deformation fields n1*n2*n3*3.
% y3     - deformation field n1*n2*n3*3.
% J1, J2 - Jacobian tensor fields n1*n2*n3*3*3.
% J3     - Jacobian tensor field n1*n2*n3*3*3.
%
% Composition of two deformations, with their Jacobian fields.
% All fields are single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT iy = spm_diffeo('invdef',y,d,M1,M2);
%
% iy - inverted deformation field of size d(1)*d(2)*d(3)*3.
% y  - original deformation field.
% M1 - An affine mapping from mm to voxels in the co-ordinate
%      system of the inverse deformation field.
% M2 - An affine mapping from voxels to mm in the co-ordinate
%      system of the forward deformation field.
%
% Inversion of a deformation field.
%
% The field is assumed to consist of a piecewise affine transformations,
% whereby each cube jointing 8 neighbouring voxels contains eight
% tetrahedra.  The mapping within each tetrahedron is assumed to be
% affine.
%
%  Reference:
%    J. Ashburner, J. Andersson and K. J. Friston (2000).
%    "Image Registration using a Symmetric Prior - in Three-Dimensions".
%    Human Brain Mapping 9(4):212-225 (appendix).
%_______________________________________________________________________
%
% FORMAT [f,dfx,dfy,dfz] = spm_diffeo('bsplins', c, y,d)
% c          - input image(s) of B-spline coefficients n1*n2*n3*n4
%              - see 'bsplinc'
% y          - points to sample n1*n2*n3*3
% d(1:3)     - degree of B-spline (from 0 to 7) along different dimensions
%              - these must be same as used by 'bsplinc'
% d(4:6)     - 1/0 to indicate wrapping along the dimensions
%
% f           - output image n1*n2*n3*n4
% dfx,dfy,dfz - sampled first derivatives
%
% c, f and y are single precision floating point.
%
% This function takes B-spline basis coefficients from spm_bsplinc,
% and re-convolves them with B-splines centred at the new sample points.
%
% Note that nearest neighbour interpolation is used instead of 0th
% degree B-splines, and the derivatives of trilinear interpolation are
% returned insted of those of 1st degree B-splines.  The difference is
% extremely subtle.
%
% c, f and y are single precision floating point.
%
%  References:
%    M. Unser, A. Aldroubi and M. Eden.
%    "B-Spline Signal Processing: Part I-Theory,"
%    IEEE Transactions on Signal Processing 41(2):821-832 (1993).
%
%    M. Unser, A. Aldroubi and M. Eden.
%    "B-Spline Signal Processing: Part II-Efficient Design and
%    Applications,"
%    IEEE Transactions on Signal Processing 41(2):834-848 (1993).
%
%    M. Unser.
%    "Splines: A Perfect Fit for Signal and Image Processing,"
%    IEEE Signal Processing Magazine, 16(6):22-38 (1999)
%
%    P. Thevenaz and T. Blu and M. Unser.
%    "Interpolation Revisited"
%    IEEE Transactions on Medical Imaging 19(7):739-758 (2000).
%
%_______________________________________________________________________
%
% FORMAT c = spm_diffeo('bsplinc',f,d)
%   f - an image
%   d(1:3) - degree of B-spline (from 0 to 7) along different dimensions
%       d(4:6) - 1/0 to indicate wrapping along the dimensions
%   c - returned volume of B-spline coefficients
%
% This function deconvolves B-splines from f, returning
% coefficients, c.  These coefficients are then passed to 'bsplins'
% in order to sample the data using B-spline interpolation.
%
%_______________________________________________________________________
%
% FORMAT f2 = spm_diffeo('samp', f1, y)
% f1 - input image(s) n1*n2*n3*n4
% y  - points to sample n1*n2*n3*3
% f2 - output image n1*n2*n3*n4
%
% Sample a function according to a deformation using trilinear interp.
% f2 = f1(y)
% f1, f2 and y are single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT f2 = spm_diffeo('push', f1, y)
% f1 - input image(s) n1*n2*n3*n4
% y  - points to sample n1*n2*n3*3
% f2 - output image n1*n2*n3*n4
%
% Push values of a function according to a deformation.  Note that the
% deformation should be the inverse of the one used with 'samp' or 'bsplins'.
% f1, f2 and y are single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT f2 = spm_diffeo('pushc', f1, y)
% f1 - input image(s) n1*n2*n3*n4
% y  - points to sample n1*n2*n3*3
% f2 - output image n1*n2*n3*n4
%
% Push values of a function according to a deformation, but using
% circulant boundary conditions.  Data wraps around.
% f1, f2 and y are single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT ut = spm_diffeo('pushg', u0, y)
% u0 - input momentum n1*n2*n3*3
% y  - points to sample n1*n2*n3*3
% ut - output momentum n1*n2*n3*3
%
% FORMAT ut = spm_diffeo('pushg', u0, y)
% u0 - input momentum n1*n2*n3*3
% y  - points to sample n1*n2*n3*3
% J  - Jacobian tensor field of y n1*n2*n3*3*3
% ut - output momentum n1*n2*n3*3
%
% Push values of a momentum field according to a deformation using
% circulant boundary conditions.  This essentially computes
% (Ad_y)^* u = |det dy| (dy)^T u(y), which is a key to the
% EPdiff equations used for geodesic shooting.
% u0, ut and y are single precision floating point.
%
%_______________________________________________________________________
%
% FORMAT f2 = spm_diffeo('resize', f1, dim)
% f1  - input fields n1*n2*n3*n4
% f2  - output field dim1*dim2*dim3*n4
% dim - output dimensions
%
% Resize a field according to dimensions dim.  This is a component of
% the multigrid approach, and is used for prolongation.
%
%_______________________________________________________________________
%
% FORMAT v2 = spm_diffeo('restrict', v1)
% v1  - input fields n1*n2*n3*n4
% v2  - output field dim1*dim2*dim3*n4
%
% Restricts a field such that its dimensions are approximately half
% their original.  This is a component of the multigrid approach.
%
%_______________________________________________________________________
%
% FORMAT J = spm_diffeo('def2jac',y)
% y - Deformation field
% J - Jacobian tensor field of y
%
% Compute Jacobian tensors from a deformation.
%
%_______________________________________________________________________
%
% FORMAT J = spm_diffeo('def2det',y)
% y - Deformation field
% j - Jacobian determinant field of y
%
% Compute Jacobian determinants from a deformation.
%
%_______________________________________________________________________
%
% FORMAT j = spm_diffeo('det',J)
% J - Jacobian tensor field
% j - Jacobian determinant field
%
% Compute determinants of Jacobian tensors.
%
%_______________________________________________________________________
%
% FORMAT dv = spm_diffeo('div',v)
% v  - velocity field
% dv - divergences of velocity field
%
% Computes divergence from velocity field.  This is indicative of rates
% of volumetric expansion/contraction.
%
%_______________________________________________________________________
%
% FORMAT [y,J] = spm_diffeo('smalldef',v,s)
% v - velocity field
% s - scaling factor
% y - small deformation
% J - approximate Jacobian tensors of small deformation (computed via
%     a matrix exponsntial of the Jacobians of the velocity field).
%
% This function is used for each time step of geodesic shooting.  It may
% change in future to use some form of Pade approximation of the
% small deformation.
%
%_______________________________________________________________________
%
% FORMAT v3 = spm_diffeo('brc', v1, v2)
% v1, v2, v3 - flow fields n1*n2*n3*3
%
% Lie Bracket.  Useful for many things
% e.g. Baker-Campbell-Haussdorf series expansion.
% The Lie bracket is denoted by
% v3 = [v1,v2]
% and on scalar fields, is computed by
% v3 = J1*v2 - J2*v1, where J1 and J2 are the Jacobian
% tensor fields. For matrices, the Lie bracket is simply
% [A,B] = A*B-B*A
%
%_______________________________________________________________________
%
% FORMAT v = spm_diffeo('dartel',v,g,f,param)
% v     - flow field n1*n2*n3*3 (single precision float)
% g     - first image n1*n2*n3*n4 (single precision float)
% f     - second image n1*n2*n3*n4 (single precision float)
% param - 9 parameters (settings)
%         - [1] Regularisation type, can take values of
%           - 0 Linear elasticity
%           - 1 Membrane energy
%           - 2 Bending energy
%         - [2][3][4] Regularisation parameters
%           - For "membrane energy", the parameters are
%             lambda, unused and id.
%           - For "linear elasticity", the parameters are
%             mu, lambda, and id
%           - For "bending energy", the parameters are
%             lambda, id1 and id2, such that regularisation is by
%             (-lambda*\grad^2 + id1)^2 + id2
%         - [5] Levenberg-Marquardt regularisation
%         - [6] Number of Full Multigrid cycles
%         - [7] Number of relaxation iterations per cycle
%         - [8] K, such that 2^K time points are used to
%               generate the deformations.  A value of zero
%               indicates a small deformation model.
%         - [9] code of 0, 1 or 2.
%               0 - asymmetric sums of squares objective function.
%               1 -  symmetric sums of squares objective function.
%               2 - assumes multinomial distribution, where template
%                   encodes the means and interpolation of template
%                   done using logs and softmax function.
%
% This is for performing a single iteration of the Dartel optimisation.
% All velocity fields and images are represented by single precision floating
% point values. Images can be scalar fields, in which case the objective
% function is the sum of squares difference.  Alternatively, images can be
% vector fields, in which case the objective function is the sum of squares
% difference between each scalar field + the sum of squares difference
% between one minus the sum of the scalar fields.
%
%_______________________________________________________________________
%
% FORMAT [y,J] = spm_diffeo('Exp', v, param)
% v - flow field
% J - Jacobian. Usually a tensor field of Jacobian matrices, but can
%     be a field of Jacobian determinants.
% param - 2 (or 3) parameters.
%         [1] K, the number of recursions (squaring steps), such
%             that exponentiation is done using an Euler-like
%             integration with 2^K time steps.
%         [2] a scaling parameter.
%         If there is a third parameter, and it is set to 1, then
%         the J will be the Jacobian determinants.
%
% A flow field is "exponentiated" to generate a deformation field
% using a scaling and squaring approach.  See the work of Arsigny
% et al, or Cleve Moler's "19 Dubious Ways" papers.
%
%_______________________________________________________________________
%
% Note that the boundary conditions are circulant throughout.
% Interpolation is trilinear, except for the resize function
% which uses a 2nd degree B-spline (without first deconvolving).
%
%_______________________________________________________________________
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_diffeo.m 4890 2012-09-03 15:19:46Z guillaume $

%-This is merely the help file for the compiled routine
error('spm_diffeo.c not compiled - see Makefile')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function a = tocmplx(a)
%TOCMPLX      Type cast to complex interval
%
%   z = tocmplx(a)
%

% written  10/16/98     S.M. Rump
% modified 11/30/98     S.M. Rump  improved speed
% modified 09/02/00     S.M. Rump  same as cintval
%                                  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 12/03/05     S.M. Rump  obsolete
%

  error('please replace "tocmplx" by "cintval"')

repo: OpenMAS