text:
function out = spm_run_fmri_est(job)
% Estimate parameters of a specified model
% SPM job execution function
% takes a harvested job data structure and call SPM functions to perform
% computations on the data.
% Input:
% job    - harvested job data structure (see matlabbatch help)
% Output:
% out    - computation results, usually a struct variable.
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

% $Id: spm_run_fmri_est.m 7354 2018-06-22 10:44:22Z guillaume $

%-Load SPM.mat file
%--------------------------------------------------------------------------
load(job.spmmat{1},'SPM');
if ~exist('SPM','var')
    error('The MAT-file does not contain an SPM variable.');
end
out.spmmat = job.spmmat;

%-Move to the directory where the SPM.mat file is
%--------------------------------------------------------------------------
original_dir = pwd;
cd(fileparts(job.spmmat{:}));

%==========================================================================
%                       R E M L   E S T I M A T I O N
%==========================================================================
if isfield(job.method,'Classical')

    %-ReML estimation of the model
    %----------------------------------------------------------------------
    SPM = spm_spm(SPM);

    %-Automatically set up contrasts for factorial designs
    %----------------------------------------------------------------------
    if isfield(SPM,'factor') && ~isempty(SPM.factor) && SPM.factor(1).levels > 1

        Ic = [];

        %-Generate contrasts
        %------------------------------------------------------------------
        cons = spm_design_contrasts(SPM);

        %-Create F-contrasts
        %------------------------------------------------------------------
        for i=1:length(cons)
            con       = cons(i).c;
            name      = cons(i).name;
            STAT      = 'F';
            [c,I]     = spm_conman('ParseCon',con,SPM.xX.xKXs,STAT);
            if all(I)
                DxCon = spm_FcUtil('Set',name,STAT,'c',c,SPM.xX.xKXs);
                if isempty(SPM.xCon)
                    SPM.xCon = DxCon;
                else
                    SPM.xCon(end+1) = DxCon;
                end
                Ic = [Ic length(SPM.xCon)];
            end
        end

        %-Create t-contrasts
        %------------------------------------------------------------------
        for i=1:length(cons)
            % Create a t-contrast for each row of each F-contrast
            % The resulting contrast image can be used in a 2nd-level analysis
            Fcon          = cons(i).c;
            STAT          = 'T';
            for r=1:size(Fcon,1)
                con       = Fcon(r,:);
                str       = cons(i).name;
                if strncmp('Interaction',str,11)
                    name  = ['Positive ',str,'_',int2str(r)];
                else
                    sp1   = find(isspace(str), 1);
                    name  = ['Positive',str(sp1:end),'_',int2str(r)];
                end

                [c,I]     = spm_conman('ParseCon',con,SPM.xX.xKXs,STAT);
                if all(I)
                    DxCon = spm_FcUtil('Set',name,STAT,'c',c,SPM.xX.xKXs);
                    if isempty(SPM.xCon)
                        SPM.xCon = DxCon;
                    else
                        SPM.xCon(end+1) = DxCon;
                    end
                    Ic = [Ic length(SPM.xCon)];
                end
            end
        end

        %-Estimate constrasts
        %------------------------------------------------------------------
        if ~isempty(Ic)
            spm('FnBanner','spm_contrasts.m');
            SPM = spm_contrasts(SPM,Ic);
        end

    end

    %-Residuals
    %----------------------------------------------------------------------
    if job.write_residuals
        VRes = spm_write_residuals(SPM,NaN);
    end

    %-Computation results
    %----------------------------------------------------------------------
    %out.spmvar = SPM;
    out.beta  = spm_file({SPM.Vbeta(:).fname}','path',SPM.swd);
    out.mask  = {fullfile(SPM.swd,SPM.VM.fname)};
    out.resms = {fullfile(SPM.swd,SPM.VResMS.fname)};
    if job.write_residuals
        out.res = spm_file({VRes(:).fname}','path',SPM.swd);
    end
    cd(original_dir);
    fprintf('Done\n');
    return
end

if job.write_residuals
    warning('Save residuals option is only implemented for classical inference.');
    job.write_residuals = false;
end

%==========================================================================
%        B A Y E S I A N   2nd   L E V E L   E S T I M A T I O N
%==========================================================================
if isfield(job.method,'Bayesian2')

    SPM = spm_spm_Bayes(SPM);

    %out.spmvar = SPM;
    cd(original_dir);
    fprintf('Done\n');
    return
end

%==========================================================================
%        B A Y E S I A N   1st   L E V E L   E S T I M A T I O N
%==========================================================================

%-Analyse specific slices or whole volume
%--------------------------------------------------------------------------
switch char(fieldnames(job.method.Bayesian.space))
  case 'volume'
      SPM.PPM.space_type  = 'volume';
      SPM.PPM.block_type  = lower(job.method.Bayesian.space.volume.block_type);
  case 'slices'
      SPM.PPM.space_type  = 'slices';
      SPM.PPM.AN_slices   = job.method.Bayesian.space.slices.numbers;
      SPM.PPM.block_type  = lower(job.method.Bayesian.space.slices.block_type);
  case 'clusters'
      SPM.PPM.space_type  = 'clusters';
      SPM.PPM.clustermask = job.method.Bayesian.space.clusters.mask;
      SPM.PPM.block_type  = lower(job.method.Bayesian.space.clusters.block_type);
  otherwise
      SPM.PPM.space_type  = 'volume';
      SPM.PPM.block_type  = 'slices';
end

%-Regression coefficient priors
%--------------------------------------------------------------------------
switch job.method.Bayesian.signal
    case 'UGL'
        SPM.PPM.priors.W  = 'Spatial - UGL';
    case 'GMRF'
        SPM.PPM.priors.W  = 'Spatial - GMRF';
    case 'LORETA'
        SPM.PPM.priors.W  = 'Spatial - LORETA';
    case 'WGL'
        SPM.PPM.priors.W  = 'Spatial - WGL';
    case 'Global'
        SPM.PPM.priors.W  = 'Voxel - Shrinkage';
    case 'Uninformative'
        SPM.PPM.priors.W  = 'Voxel - Uninformative';
    otherwise
        error('Unknown prior for W in Bayesian 1st level estimation.');
end

%-Number of AR coefficients
%--------------------------------------------------------------------------
SPM.PPM.AR_P              = job.method.Bayesian.ARP;

%-AR coefficient priors
%--------------------------------------------------------------------------
if isfield(job.method.Bayesian.noise,'UGL')
    SPM.PPM.priors.A      = 'Spatial - UGL';
elseif isfield(job.method.Bayesian.noise,'GMRF')
    SPM.PPM.priors.A      = 'Spatial - GMRF';
elseif isfield(job.method.Bayesian.noise,'LORETA')
    SPM.PPM.priors.A      = 'Spatial - LORETA';
elseif isfield(job.method.Bayesian.noise,'tissue_type')
    SPM.PPM.priors.A      = 'Discrete';
    SPM.PPM.priors.SY     = char(job.method.Bayesian.noise.tissue_type);
elseif isfield(job.method.Bayesian.noise,'Robust')
    SPM.PPM.priors.A      = 'Robust';
    SPM.PPM.AR_P          = 0;
    SPM.PPM.update_F      = 1;
end

%-Define an empty contrast
%--------------------------------------------------------------------------
NullCon      = spm_FcUtil('Set','','P','c',[],1);
NullCon.X0   = [];
NullCon.iX0  = [];
NullCon.X1o  = [];
NullCon.eidf = 1;

SPM.xCon     = [];
%-Set up contrasts for 2nd-level ANOVA
%--------------------------------------------------------------------------
if strcmp(job.method.Bayesian.anova.second,'Yes') && isfield(SPM,'factor')
    cons = spm_design_contrasts(SPM);
    for i=1:length(cons)
        % Create a simple contrast for each row of each F-contrast
        % The resulting contrast image can be used in a 2nd-level analysis
        Fcon = cons(i).c;
        for r=1:size(Fcon,1)
            % Normalise contrast st. sum of positive elements is 1
            % and sum of negative elements is 1
            con = Fcon(r,:);
            con = con / length(find(con==1));

            % Change name
            str = cons(i).name;
            sp1 = find(str==' ', 1);
            if strcmp(str(1:11),'Interaction')
                name = ['Positive ',str,'_',int2str(r)];
            else
                name = ['Positive',str(sp1:end),'_',int2str(r)];
            end

            DxCon = NullCon;
            DxCon.name = name;
            DxCon.c = con';

            if isempty(SPM.xCon)
                SPM.xCon = DxCon;
            else
                SPM.xCon(end+1) = DxCon;
            end
        end
    end
end

%-Set up user-specified simple contrasts
%--------------------------------------------------------------------------
K               = size(SPM.xX.X,2);
for c = 1:length(job.method.Bayesian.gcon)
    DxCon       = NullCon;
    DxCon.name  = job.method.Bayesian.gcon(c).name;
    convec      = job.method.Bayesian.gcon(c).convec(:);
    if length(convec) == K
        DxCon.c = convec;
    elseif length(convec) < K
        fprintf('Zero padding.');                                       %-#
        DxCon.c = [convec; zeros(K-length(convec),1)];
    else
        warning(['User-specified contrast nb %d has %d entries '...
            'but there are %d regressors - ignored.'], c,length(convec),K);
        DxCon   = [];
    end

    if isempty(SPM.xCon)
        SPM.xCon = DxCon;
    else
        SPM.xCon(end+1) = DxCon;
    end
end

%-Compute log evidence maps
%--------------------------------------------------------------------------
if strcmp(job.method.Bayesian.LogEv,'Yes')
    SPM.PPM.update_F      = 1;
    SPM.PPM.compute_det_D = 1;
end

%-1st level Bayesian ANOVA ?
%--------------------------------------------------------------------------
bayes_anova = 0;
if strcmp(job.method.Bayesian.anova.first,'Yes')
    bayes_anova           = 1;
    SPM.PPM.update_F      = 1; % Compute evidence for each model
    SPM.PPM.compute_det_D = 1;
end

%-Variational Bayes estimation
%--------------------------------------------------------------------------
SPM = spm_spm_vb(SPM);

%-Bayesian ANOVA using model comparison
%--------------------------------------------------------------------------
if bayes_anova
    % We don't want to estimate contrasts for each different model
    SPM.xCon = [];
    spm_vb_ppm_anova(SPM);
end

%out.spmvar = SPM;
cd(original_dir);
fprintf('Done\n')

author: spm


text:
function [K, sK, n1] = dexpKernCompute(kern, x1, x2)

% DEXPKERNCOMPUTE Compute the double exponential kernel,
%
% k(x_i, x_j) = 0.5 * sigma2 * theta * exp(-theta*abs(x_i - x_j)),
%
% given the parameters (theta and sigma), and t.
%
% FORMAT
% DESC computes the kernel parameters for the double exponential kernel
% given the input matrices associated with rows and columns.
% ARG kern : the kernel structure for which the kernel matrix is computed.
% ARG x1 : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN K : the kernel matrix computed at the given points.
% RETURN sK: normalised kernel matrix (i.e. variance set to 1).
% RETURN n1: L1 distance between each row of x1 and x2.
%
% FORMAT
% DESC computes the kernel matrix for the double exponential kernel
% given a matrix of inputs.
% ARG kern : the kernel structure for which the kernel matrix is computed.
% ARG x1 : the input matrix associated with the rows and the columns.
% RETURN K : the kernel matrix computed at the given points.
% RETURN sK: normalised kernel matrix (i.e. no multiplication by the decay
% or variance of the exponential).
% RETURN n1: L1 distance between each row of x1 and x2.
%
% SEEALSO : ouKernParamInit, kernCompute, kernCreate, ouKernDiagCompute,
%
% COPYRIGHT : David Luengo, 2009
%
% COPYRIGHT : Neil D. Lawrence, 2009

% KERN

if nargin < 3
  n1 = sqrt(dist2(x1, x1));
else
  n1 = sqrt(dist2(x1, x2));
end

sK = 0.5 * exp(-kern.decay * n1);
K = sK * kern.decay * kern.variance;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [tUs,odrIdx,TXmean,Wgt,LDAU] = MPCALDA(TX,gndTX,MPCADADim,testQ,maxK)
% MPCA+LDA: Multilinear Principle Component Analysis plus Linear
% Discriminant Analysis
%
% %[Prototype]%
% function [tUs,odrIdx,TXmean,Wgt,LDAU]=MPCALDA(TX,gndTX,MPCADADim,testQ,maxK)
%
% %[Author Notes]%
% Author: Haiping LU
% Email : hplu@ieee.org   or   eehplu@gmail.com
% Affiliation: Department of Electrical and Computer Engineering
%              University of Toronto
% Release date: June 24, 2008 (Version 1.1)
% Please email me if you have any problem, question, or suggestion
%
% %[Algorithm]%:
% This function implements the Multilinear Principal Component Analysis
% plus Linear Discriminant Analysis (MPCA+LDA) algorithm presented in the
% follwing paper:
%    Haiping Lu, K.N. Plataniotis, and A.N. Venetsanopoulos,
%    "MPCA: Multilinear Principal Component Analysis of Tensor Objects",
%    IEEE Transactions on Neural Networks,
%    Vol. 19, No. 1, Page: 18-39, January 2008.
% Please reference this paper when reporting work done using this code.
%
% %[Toolbox needed]%:
% This function needs the tensor toolbox available at
% http://csmr.ca.sandia.gov/~tgkolda/TensorToolbox/
%
% %[Syntax]%:
%    [tUs,odrIdx,TXmean,Wgt,LDAU] = MPCALDA(TX,gndTX,MPCADADim,testQ,maxK)
%
% %[Inputs]%:
%    TX: the input training data in tensorial representation, the last mode
%        is the sample mode. For Nth-order tensor data, TX is of
%        (N+1)th-order with the (N+1)-mode to be the sample mode.
%        E.g., 30x20x10x100 for 100 samples of size 30x20x10
%        If your training data is too big, resulting in the "out of memory"
%        error, you could work around this problem by reading samples one
%        by one from the harddisk, or you could email me for help.
%
%    gndTX: the ground truth class labels (1,2,3,...) for the training data
%           E.g., a 100x1 vector if there are 100 samples
%
%    MPCADADim: the number of MPCA features to be fed into LDA.
%               This parameter could have a significant effect on the
%               recognition performance, as in PCA+LDA. It is strongly
%               suggested that you should test different numbers of
%               MPCADADim to get the best performance.
%
%    testQ: the percentage of variation kept in each mode, suggested value
%           is 97, and you can try other values, e.g., from 95 to 100, to
%           see whether better performance can be obtained.
%
%    maxK: the maximum number of iterations, suggested value is 1, and you
%          can try a larger value if computational time is not a concern.
%
% %[Outputs]%:
%    tUs: the multilinear projection, consiting of N
%         projection matrices, one for each mode
%
%    odrIdx: the ordering index of projected features in decreasing
%            discriminality for vectorizing the projected tensorial
%            features
%
%    TXmean: the mean of the input training samples TX
%
%    Wgt: the weight vector for use in modified distance measures. Please
%         refer to Section IV.C of the paper.
%
%    LDAU: the standard LDA projection matrix
%
% %[Supported tensor order]%
% This function supports N=2,3,4, for other order N, please modify the
% codes accordingly or email hplu@ieee.org or eehplu@gmail.com for help
%
% %[Examples]%
%%%%%%%%%%%%%%%%%%%%%%%%%%Example on 2D face data%%%%%%%%%%%%%%%%%%%%%%%%%%
%       load FERETC80A45%each sample is a second-order tensor of size 32x32
%       N=ndims(fea2D)-1;%Order of the tensor sample
%       Is=size(fea2D);%32x32x320
%       numSpl=Is(3);%There are 320 face samples
%       MPCADADim=80;%80 most discriminative MPCA features for LDA
%       testQ=97;%Keep 97% variation in each mode
%       maxK=1;%One iteration only
%       [tUs,odrIdx,TXmean,Wgt,LDAU] = MPCALDA(fea2D,gnd,MPCADADim,testQ,maxK);
%       fea2Dctr=fea2D-repmat(TXmean,[ones(1,N), numSpl]);%Centering
%       newfea=ttm(tensor(fea2Dctr),tUs,1:N);%MPCA projection
%       %Vectorization of the tensorial feature
%       newfeaDim=size(newfea,1)*size(newfea,2);
%       newfea=reshape(newfea.data,newfeaDim,numSpl)';%Note: Transposed
%       selfea = newfea(:,odrIdx);%Take only the first MPCADADim ones.
%       %Note that the length of odrIdx is already "MPCADADim"
%       LDAfea=selfea*LDAU;%LDA projection
%       %LDAfea is the MPCA+LDA feature vector for input to classical
%       %classifiers (e.g., nearest neighbor classifier).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%Example on 3D gait data%%%%%%%%%%%%%%%%%%%%%%%%%%
%       load USF17Gal %each sample is a third-order tensor of size 32x22x10
%       N=ndims(fea3D)-1;%Order of the tensor sample
%       Is=size(fea3D);%32x22x10x731
%       numSpl=Is(4);%There are 731 gait samples
%       MPCADADim=200;%200 most discriminative MPCA features for LDA
%       testQ=97;%Keep 97% variation in each mode
%       maxK=1;%One iteration only
%       [tUs,odrIdx,TXmean,Wgt,LDAU] = MPCALDA(fea3D,gnd,MPCADADim,testQ,maxK);
%       fea3Dctr=fea3D-repmat(TXmean,[ones(1,N), numSpl]);%Centering
%       newfea=ttm(tensor(fea3Dctr),tUs,1:N);%MPCA projection
%       %Vectorization of the tensorial feature
%       newfeaDim=size(newfea,1)*size(newfea,2)*size(newfea,3);
%       newfea=reshape(newfea.data,newfeaDim,numSpl)';%Note: Transposed
%       selfea = newfea(:,odrIdx);%Take only the first MPCADADim ones.
%       %Note that the length of odrIdx is already "MPCADADim"
%       LDAfea=selfea*LDAU;%LDA projection
%       %LDAfea is the MPCA+LDA feature vector for input to classical
%       %classifiers (e.g., nearest neighbor classifier).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %[Notes]%:
% A. Developed using Matlab R2006a
% B. Revision history:
%       Version 1.0 released on March 1, 2008
%       Version 1.1 released on June 24, 2008
%           ---Example usage on 2D data is included
%           ---Recommendation on parameter MPCADADim is updated
%       Version 1.2 released on March 12, 2011
%           ---Insertion of line 275
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%TX: (N+1)-dimensional tensor of Tensor Sample Dimension x NumSamples
N=ndims(TX)-1;%The order of samples.
IsTX=size(TX);
Is=IsTX(1:N);%The dimensions of the tensor
numSpl=IsTX(N+1);%Number of samples

%%%%%%%%%%%%%Zero-Mean%%%%%%%%%%
TXmean=mean(TX,N+1);%The mean
TX=TX-repmat(TXmean,[ones(1,N), numSpl]);%Centering

%The full projection for initialization
Qs=ones(N,1)*testQ;
Us=cell(N,1);
tUs=cell(N,1);
Lmds=cell(N,1);
for n=1:N
    In=Is(n);Phi=zeros(In,In);
    for m=1:numSpl
        switch N
            case 2
                Xm=TX(:,:,m);
            case 3
                Xm=TX(:,:,:,m);
            case 4
                Xm=TX(:,:,:,:,m);
            otherwise
                error('Order N not supported. Please modify the code here or email hplu@ieee.org for help.')
        end
        tX=tensor(Xm);
        tXn=tenmat(tX,n);
        Xn=tXn.data;
        Phi=Phi+Xn*Xn';
    end
    [Un,Lmdn]=eig(Phi);
    Lmd=diag(Lmdn);
    [stLmd,stIdx]=sort(Lmd,'descend');
    Us{n}=Un(:,stIdx);
    tUs{n}=Us{n}';
    Lmds{n}=Lmd(stIdx);
end

%Cumulative distribution of eigenvalues
cums=cell(N,1);
for n=1:N
    In=length(Lmds{n});
    cumLmds=zeros(In,1);
    Lmd=Lmds{n};
    cumLmds(1)=Lmd(1);
    for in=2:In
        cumLmds(in)=cumLmds(in-1)+Lmd(in);
    end
    cumLmds=cumLmds./sum(Lmd);
    cums{n}=cumLmds;
end

%MPCA Iterations
if maxK>0
    tPs=cell(N,1);
    pUs=cell(N,1);
    %%%%%%%%%%%%%Determine Rn, the dimension of projected space%%%%
    for n=1:N
        cum=cums{n};
        idxs=find(cum>=Qs(n)/100);
        Ps(n)=idxs(1);
        tUn=tUs{n};
        tPn=tUn(1:Ps(n),:);
        tPs{n}=tPn;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for iK=1:maxK
        for n=1:N
            In=Is(n);
            Phi=double(zeros(In,In));
            for m=1:numSpl
                switch N
                    case 2
                        Xm=TX(:,:,m);
                    case 3
                        Xm=TX(:,:,:,m);
                    case 4
                        Xm=TX(:,:,:,:,m);
                    otherwise
                        error('Order N not supported. Please modify the code here or email hplu@ieee.org for help.')
                end
                tX=ttm(tensor(Xm),tPs,-n);
                tXn=tenmat(tX,n);
                Xn=tXn.data;
                Phi=Phi+Xn*Xn';
            end
            Pn=Ps(n);
            Phi=double(Phi);
            if Pn<In
                option=struct('disp',0);
                [pUs{n},pLmdn]=eigs(Phi,Pn,'lm',option);
                pLmds{n}=diag(pLmdn);
            else
                [pUn,pLmdn]=eig(Phi);
                pLmd=diag(pLmdn);
                [stLmd,stIdx]=sort(pLmd,'descend');
                pUs{n}=pUn(:,stIdx(1:Pn));
                pLmds{n}=pLmd(stIdx(1:Pn));
            end
            tPs{n}=pUs{n}';
        end
    end
    Us=pUs;
    tUs=tPs;
    Lmds=pLmds;
    Is=Ps;
else
    if testQ<100
        error('At least one iteration is needed');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Yps=ttm(tensor(TX),tUs,1:N);%MPCA projections of samples TX
vecDim=1;
for n=1:N, vecDim=vecDim*Is(n); end
vecYps=reshape(Yps.data,vecDim,numSpl); %vectorization of Yps
%%%%%%%%%%%%%%Now vecYps contains the feature vectors for training data

%%%%%%%%%%%%%%%Sort according to Fisher's discriminality%%%%%%%%%%%%%%%
classLabel = unique(gndTX);
nClass = length(classLabel);%Number of classes
ClsIdxs=cell(nClass);
Ns=zeros(nClass,1);
for i=1:nClass
    ClsIdxs{i}=find(gndTX==classLabel(i));
    Ns(i)=length(ClsIdxs{i});
end
Ymean=mean(vecYps,2);
TSW=zeros(vecDim,1);
TSB=zeros(vecDim,1);
for i=1:nClass
    clsYp=vecYps(:,ClsIdxs{i});
    clsMean=mean(clsYp,2);
    FtrDiff=clsYp-repmat(clsMean,1,Ns(i));
    TSW=TSW+sum(FtrDiff.*FtrDiff,2);
    meanDiff=clsMean-Ymean;
    TSB=TSB+Ns(i)*meanDiff.*meanDiff;
end
FisherRatio=TSB./TSW;
[stRatio,odrIdx]=sort(FisherRatio,'descend');

%%%%%%%%%%%%%%%%%%%%%%%%%%LDA%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if length(odrIdx)<MPCADADim, MPCADADim=length(odrIdx);end
odrIdx=odrIdx(1:MPCADADim);%Take only the first a few
vecYps=vecYps(odrIdx,:);%Input to LDA
Ymean = mean(vecYps,2);%should be  zero
SB = zeros(MPCADADim);
for i = 1:nClass
    clsYp=vecYps(:,ClsIdxs{i});
    clsMean=mean(clsYp,2)-Ymean;
    SB = SB + Ns(i)*clsMean*clsMean';
end
SW = vecYps*vecYps' - SB;
if rank(SW)<MPCADADim
    SW=SW+1e-6*eye(MPCADADim);
end
option=struct('disp',0);
[LDAU, LDAV] = eigs(inv(SW)*SB,nClass-1,'lm',option);
%Calculate the weight
LDAV=diag(LDAV);
Wgt=sqrt(LDAV);

author: Sable

repo: mcbench-benchmarks


text:
function out = feval(F, x, varargin)
%FEVAL   Evaluate a CHEBFUN.
%   FEVAL(F, X) evaluates a CHEBFUN F at the points in X.  If F is a quasimatrix
%   with columns F1, ..., FN, then the result will be [F1(X), ..., FN(X)], the
%   horizontal concatenation of the results of evaluating each column at the
%   points in X.
%
%   FEVAL(F, 'left'), FEVAL(F, 'start'), and FEVAL(F, '-') return the value of F
%   at the left endpoint of its domain.  FEVAL(F, 'right'), FEVAL(F, 'end'), and
%   FEVAL(F, '+') do the same for the right endpoint.
%
%   FEVAL(F, X, 'left') and FEVAL(F, X, '-') evaluate F at the points in X,
%   using left-hand limits to evaluate F at any breakpoints. FEVAL(F, X,
%   'right') and FEVAL(F, X, '+') do the same but using right-hand limits.
%
%   F(X), F('left'), F(X, 'left'), etc, are equivalent syntaxes.
%
%   Example:
%     f = chebfun(@(x) 1./(1 + 25*x.^2));
%     y = feval(f, linspace(-1, 1, 100).');
%
% See also SUBSREF.

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

% If F or x are empty, there's nothing to do.
if ( isempty(F) )
    out = [];
    return
elseif ( isempty(x) )
    % Return empty matrix with dimensions of the appropriate size.
    out = zeros(size(x));
    return
end

% Deal with FEVAL(FUNCTION_HANDLE, CHEBFUN) type calls.
if ( isa(F, 'function_handle') )
    out = F(x, varargin{:});
    return
end

% Deal with FEVAL(CHEBFUN, CHEBFUN) type calls.
if ( isa(x, 'chebfun') )
    out = compose(x, F);
    return
end

%% LEFT / RIGHT / END VALUES:
% Support for feval(f, 'left') and feval(f, 'end'), etc.
if ( ischar(x) )
    if ( any(strcmpi(x, {'left', 'start' ,'-'})) )
        out = F.pointValues(1,:);
    elseif ( any(strcmpi(x, {'right', 'end', '+'})) )
        out = F.pointValues(end,:);
    else
        error('CHEBFUN:CHEBFUN:feval:strInput', ...
            'Unknown input argument "%s".', x);
    end
    if ( F(1).isTransposed )
        out = out.';
    end
    return
end

%% EVALUATE EACH COLUMN:
if ( numel(F) == 1 )
    out = columnFeval(F, x, varargin{:});
else
    % Deal with quasimatrices:
    out = cell(1, numel(F));
    for k = 1:numel(F)
        out{k} = columnFeval(F(k), x, varargin{:});
    end
    out = cell2mat(out);
end

%% DEAL WITH ROW CHEBFUNS:
if ( F(1).isTransposed )
    % We got a passed a row CHEBFUN. If X had more than two dimensions, we can't
    % simply transpose the output from above, instead, we need to use permute.
    ndimsx = ndims(x);
    if ( ndimsx <= 2 )
        out = out.';
    else
        % We define "transposition" in this case to mean the switching of the
        % first two dimensions.
        out = permute(out, [2 1 3:ndimsx]);
    end
end

end

function out = columnFeval(f, x, varargin)
% Act on each column. Note that columnFeval ignores the transpose state of F.

%% INITIALISE:

% Reshape x to be a column vector.
sizex = size(x);
ndimsx = ndims(x);
x = x(:);

% Initialise output:
numFuns = numel(f.funs);

funs = f.funs;
dom = f.domain;
isTrans = f.isTransposed;

%% LEFT AND RIGHT LIMITS:
% Deal with feval(f, x, 'left') and feval(f, x, 'right'):
leftFlag = 0; rightFlag = 0;
if ( nargin > 2 )
    lr = varargin{1};
    leftFlag = any(strcmpi(lr, {'left', '-'}));
    rightFlag = any(strcmpi(lr, {'right', '+'}));
    if ( ~(leftFlag || rightFlag) )
        if ( ischar(lr) )
            error('CHEBFUN:CHEBFUN:feval:leftRightChar',...
                'Unknown input argument "%s".', lr);
        else
            error('CHEBFUN:CHEBFUN:feval:leftRight', 'Unknown input argument.');
        end
    end
end

%% VALUES FROM FUNS:

if ( numFuns == 1 )

    % Things are simple when there is only a single FUN:
    out = feval(funs{1}, x(:), varargin{:});
    numCols = size(out, 2);

else

    % For multiple FUNs we must determine which FUN corresponds to each x.

    % Initialise output matrix:
    numCols = numColumns(f);
    out = zeros(numel(x), numCols);

    % Replace the first and last domain entries with +/-inf. (Since we want to
    % use FUN{1} if real(x) < dom(1) and FUN{end} if real(x) > dom(end)).
    domInf = [-inf, dom(2:end-1), inf];

    % Loop over each FUN. If real(x) is in [dom(k) dom(k+1)] then use FUN{k}.
    xReal = real(x);
    for k = 1:numFuns
        I = ( xReal >= domInf(k) ) & ( xReal < domInf(k+1) );
        if ( any(I(:)) )
            % Evaluate the appropriate fun:
            out(I,:) = feval(funs{k}, x(I), varargin{:});
        end
    end

end

%% POINTVALUES:
% If the evaluation point corresponds to a breakpoint, we get the value from

% f.pointValues. However, if the 'left' or 'right' flag is given, we first
% modify the entry in point values to take either the left- or right-sided
% limit, respectively.

[xAtBreaks, domIndices] = ismember(x, dom);
if ( any(xAtBreaks) )
    if ( leftFlag )
        % Note that for leftFlag we use 'rval-local', which corresponds to the
        % function value at the right part of the subdomain.
        rvals = get(f, 'rval-local'); if ( isTrans ), rvals = rvals.'; end
        pointValues = [f.pointValues(1,:); rvals];
    elseif ( rightFlag )
        % Similarly rightFlag uses lval-local.
        lvals = get(f, 'lval-local'); if ( isTrans ), lvals = lvals.'; end
        pointValues = [lvals; f.pointValues(end,:)];
    else
        pointValues = f.pointValues;
    end
    % Set the output values for any values of x at the breakpoints.
    out(xAtBreaks,:) = pointValues(domIndices(xAtBreaks),:);
end

%% RESHAPE FOR OUTPUT:
% Reshape fx, which is a column vector or horizontal concatenation of column
% vectors, to be of the appropriate size, and handle transposition.
sizefx = sizex;
sizefx(2) = numCols*sizex(2);
if ( ndimsx == 2 )
    % If x was just a matrix or vector, the reshape is straightforward.
    out = reshape(out, sizefx);
else
    % If x had more than two dimensions, we have to be more careful.  The
    % cell2mat(mat2cell(...).') effects a transpose which keeps certain
    % columnar blocks of the fx matrix intact, putting the entries in the
    % correct order for reshape().
    blockLength = sizex(1)*sizex(2);
    blocksPerCol = prod(sizex(3:end));
    out = reshape(cell2mat(mat2cell(out, blockLength*ones(1, blocksPerCol), ...
        ones(1, numCols)).'), sizefx);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [y]=multifuncrs(X, funs, eps, varargin)
% Cross approximation of a (vector-)function of several TT-tensors.
%   [Y]=MULTIFUNCRS(X,FUNS,EPS, VARARGIN)
%   Computes approximation to the functions FUNS(X{1},...,X{N}) with accuracy EPS
%   X should be a cell array of nx TT-tensors of equal sizes.
%   The function FUNS should receive a 2d array V of sizes I x N, where the
%   first dimension stays for the reduced set of spatial indices, and  the
%   second is the enumerator of X.
%   The returned sizes should be I x D2, where D2 is the number of
%   components in FUNS. D2 should be either provided as the last (d+1)-th
%   TT-rank of the initial guess, or given explicitly as an option (see
%   below).
%   For example, a linear combination reads FUNS=@(x)(x*W), W is a N x D2
%   matrix.
%
%   Options are provided in form
%   'PropertyName1',PropertyValue1,'PropertyName2',PropertyValue2 and so
%   on. The parameters are set to default (in brackets in the following)
%   The list of option names and default values are:
%       o y0 - initial approximation [random rank-2 tensor]
%       o nswp - maximal number of DMRG sweeps [10]
%       o rmax - maximal TT rank [Inf]
%       o verb - verbosity level, 0-silent, 1-sweep info, 2-block info [1]
%       o kickrank - the rank-increasing parameter [5]
%       o d2 - the last rank of y, that is dim(FUNS) [1]
%       o qr - do (or not) qr before maxvol [false]
%
%   The method is based on the alternating approximation, with
%   the one-block enrichment via KICKRANK random vectors or randomized AMR.
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

nswp = 10;
kickrank = 5;
y = [];
verb = 1;
% kicktype = 'rand';
kicktype = 'amr-two';
pcatype = 'svd';
%pcatype = 'uchol';
rmax = Inf;
d2 = 1;
wasrand = false;
trunctype = 'fro';
do_qr = false;
% trunctype = 'cross';

for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'nswp'
            nswp=varargin{i+1};
        case 'y0'
            y=varargin{i+1};
        case 'kickrank'
            kickrank=varargin{i+1};
        case 'rmax'
            rmax=varargin{i+1};
        case 'verb'
            verb=varargin{i+1};
        case 'kicktype'
            kicktype=varargin{i+1};
        case 'pcatype'
            pcatype=varargin{i+1};
        case 'trunctype'
            trunctype=varargin{i+1};
        case 'd2'
            d2=varargin{i+1};
        case 'qr'
            do_qr = varargin{i+1};
        otherwise
            error('Unrecognized option: %s\n',varargin{i});
    end
end

nx = numel(X);
d = X{1}.d;
n = X{1}.n;
rx = zeros(d+1,nx);
crX = cell(d,nx);
for i=1:nx
    rx(:,i) = X{i}.r;
    crX(:,i) = core2cell(X{i});
end;

if (isempty(y))
    ry = d2*ones(d+1,1); ry(1)=1;
    y = tt_rand(n, d, ry);
    wasrand = true;
end;
ry = y.r;
cry = core2cell(y);

Ry = cell(d+1,1);
Ry{1} = 1; Ry{d+1}=1;
Rx = cell(d+1,nx);
for i=1:nx
    Rx{1,i}=1; Rx{d+1,i}=1;
end;

block_order = [+(d), -(d)];

% Orth
for i=1:d-1
    cr = cry{i}; % r1,n,d2,r2
    cr = reshape(cr, ry(i)*n(i), ry(i+1));
    [cr, rv]=qr(cr, 0);
    cr2 = cry{i+1};
    cr2 = reshape(cr2, ry(i+1), n(i+1)*ry(i+2));
    cr2 = rv*cr2;
    ry(i+1) = size(cr, 2);
    cr = reshape(cr, ry(i), n(i), ry(i+1));
    cry{i+1} = reshape(cr2, ry(i+1), n(i+1), ry(i+2));
    cry{i} = cr;

    % Interface matrix for Y
    Ry{i+1} = Ry{i}*reshape(cr, ry(i), n(i)*ry(i+1));
    Ry{i+1} = reshape(Ry{i+1}, ry(i)*n(i), ry(i+1));
    if (wasrand)
        curind = [];
        while numel(curind)<ry(i+1)
            curind = [curind; ceil(rand(ry(i+1), 1)*(n(i)*ry(i)))];
            curind = unique(curind);
        end;
        curind = curind(1:ry(i+1));
    else

        curind = maxvol2(Ry{i+1},'qr',do_qr);
    end;
    Ry{i+1} = Ry{i+1}(curind, :);
    % Interface matrices for X
    for j=1:nx
        Rx{i+1,j} = reshape(crX{i,j}, rx(i,j), n(i)*rx(i+1,j));
        Rx{i+1,j} = Rx{i,j}*Rx{i+1,j};
        Rx{i+1,j} = reshape(Rx{i+1,j}, ry(i)*n(i), rx(i+1,j));
        Rx{i+1,j} = Rx{i+1,j}(curind, :);
    end;
end;

d2 = ry(d+1);
ry(d+1) = 1;
cry{d} = permute(cry{d}, [3,1,2]); % d2, rd, nd

last_sweep = false;
swp = 1;

% dy_old = ones(d,1);
dy = zeros(d,1);
max_dy = 0;
% For extra-rank addition
% dpows = ones(d,1)*min_dpow;
% dranks = zeros(d,1);

cur_order = block_order;
order_index = 2;
i = d;
dir = sign(cur_order(order_index));

% DMRG sweeps
while (swp<=nswp)||(dir>0)

    oldy = reshape(cry{i}, d2*ry(i)*n(i)*ry(i+1), 1);

    if (~last_sweep)
        % Compute the X superblocks
        curbl = zeros(ry(i)*n(i)*ry(i+1), nx);
        for j=1:nx
            cr = reshape(crX{i,j}, rx(i,j), n(i)*rx(i+1,j));
            cr = Rx{i,j}*cr;
            cr = reshape(cr, ry(i)*n(i), rx(i+1,j));
            cr = cr*Rx{i+1,j};
            curbl(:,j) = cr(:);
        end;
        % Call the function
        newy = funs(curbl); % sizes: rnr x nx -> rnr x d2
        % Multiply with inverted Ry
        newy = reshape(newy, ry(i), n(i)*ry(i+1)*d2);
        newy = (Ry{i}) \ newy;
        newy = reshape(newy, ry(i)*n(i)*ry(i+1), d2);
        newy = reshape(newy.', d2*ry(i)*n(i), ry(i+1));
        newy = newy / (Ry{i+1});
        newy = reshape(newy, d2*ry(i)*n(i)*ry(i+1), 1);
    else
        newy = oldy;
    end;

    dy(i) = norm(newy(:)-oldy)/norm(newy(:));
    max_dy = max(max_dy, dy(i));

    % Truncation
    if (dir>0) % left-to-right
        newy = reshape(newy, d2, ry(i)*n(i)*ry(i+1));
        newy = reshape(newy.', ry(i)*n(i), ry(i+1)*d2);
    else
        newy = reshape(newy, d2*ry(i), n(i)*ry(i+1));
    end;

    if (kickrank>=0)
        [u,s,v]=svd(newy, 'econ');
        s = diag(s);
        if (strcmp(trunctype, 'fro'))||(last_sweep)
            r = my_chop2(s, eps/sqrt(d)*norm(s));
        else
            % Truncate taking into account the (r+1) overhead in the cross
            cums = (s.*(2:numel(s)+1)').^2;
            cums = cumsum(cums(end:-1:1));
            cums = cums(end:-1:1)./cums(1);
            r = find(cums<(eps^2/d), 1);
            if (isempty(r))
                r = numel(s);
            end;
        end;
        r = min(r, rmax);
        r = min(r, numel(s));
    else
        if (dir>0)
            [u,v]=qr(newy, 0);
            v=v';
            r = size(u,2);
            s = ones(r,1);
        else
            [v,u]=qr(newy.', 0);
            v=conj(v);
            u=u.';
            r = size(u,2);
            s = ones(r,1);
        end;
    end;

    if (verb>1)
    	fprintf('=multifuncrs=   block %d{%d}, dy: %3.3e, r: %d\n', i, dir, dy(i), r);
    end;

    % Kicks and interfaces
    if (dir>0)&&(i<d) % left-to-right, kickrank, etc
        u = u(:,1:r);
        v = v(:,1:r)*diag(s(1:r));

        % kick
        radd = 0; rv = 1;
        if (~last_sweep)&&(kickrank>0)
            if (strcmp(kicktype, 'amr-two'))
                % AMR(two)-like kick. See also the M.Sc.Thesis by D. Zheltkov.
                % The left indices are nested, but the right are chosen
                % randomly. In Zheltkov's work, from all possible n^(d-k)
                % values. However, in the functional-cross it would result
                % in a d^2 complexity. Here, I use only the neighbouring
                % core for randomization. Additionally, the actual kick is
                % performed via the Z=PCA(supercore), since otherwise the
                % rank grows too high.

                % Compute the X superblocks
                ind2 = unique(ceil(rand(ry(i+1), 1)*(ry(i+2)*n(i+1))));
                rkick = numel(ind2);
                curbl = zeros(ry(i)*n(i)*rkick, nx);
                for j=1:nx
                    cr1 = reshape(crX{i,j}, rx(i,j), n(i)*rx(i+1,j));
                    cr1 = Rx{i,j}*cr1;
                    cr1 = reshape(cr1, ry(i)*n(i), rx(i+1,j));
                    cr2 = reshape(crX{i+1,j}, rx(i+1,j)*n(i+1), rx(i+2,j));
                    cr2 = cr2*Rx{i+2,j}; % now its size rx
                    cr2 = reshape(cr2, rx(i+1,j), n(i+1)*ry(i+2));
                    cr2 = cr2(:, ind2);
                    curbl(:,j) = reshape(cr1*cr2, ry(i)*n(i)*rkick, 1);
                end;
                % Call the function
                uk = funs(curbl); % rnr, d2
                uk = reshape(uk, ry(i), n(i)*rkick*d2);
                uk = Ry{i} \ uk;
                uk = reshape(uk, ry(i)*n(i), rkick*d2);
                if (strcmp(pcatype, 'svd'))
                    [uk,sk,vk]=svd(uk, 'econ');
                    uk = uk(:,1:min(kickrank, size(uk,2)));
                else
                    uk = uchol(uk.', kickrank+1);
                    uk = uk(:,end:-1:max(end-kickrank+1,1));
                end;
            else
                uk = rand(ry(i)*n(i), kickrank);
            end;
            [u,rv]=qr([u,uk], 0);
            radd = size(uk,2);
        end;
        v = [v, zeros(ry(i+1)*d2, radd)];
        v = rv*(v');
        r = size(u,2);

        cr2 = cry{i+1};
        cr2 = reshape(cr2, ry(i+1), n(i+1)*ry(i+2));
        v = reshape(v, r*ry(i+1), d2);
        v = reshape(v.', d2*r, ry(i+1));
        v = v*cr2; % size r+radd, n2, r3

        ry(i+1) = r;

        u = reshape(u, ry(i), n(i), r);
        v = reshape(v, d2, r, n(i+1), ry(i+2));

        % Stuff back
        cry{i} = u;
        cry{i+1} = v;

        % Recompute left interface matrices
        % Interface matrix for Y
        Ry{i+1} = Ry{i}*reshape(u, ry(i), n(i)*ry(i+1));
        Ry{i+1} = reshape(Ry{i+1}, ry(i)*n(i), ry(i+1));
        curind = maxvol2(Ry{i+1},'qr',do_qr);
        Ry{i+1} = Ry{i+1}(curind, :);
        % Interface matrices for X
        for j=1:nx
            Rx{i+1,j} = reshape(crX{i,j}, rx(i,j), n(i)*rx(i+1,j));
            Rx{i+1,j} = Rx{i,j}*Rx{i+1,j};
            Rx{i+1,j} = reshape(Rx{i+1,j}, ry(i)*n(i), rx(i+1,j));
            Rx{i+1,j} = Rx{i+1,j}(curind, :);
        end;
    elseif (dir<0)&&(i>1) % right-to-left
        u = u(:,1:r)*diag(s(1:r));
        v = conj(v(:,1:r));
        % kick
        radd = 0; rv = 1;
        if (~last_sweep)&&(kickrank>0)
            if (strcmp(kicktype, 'amr-two'))
                % Compute the X superblocks
                ind2 = unique(ceil(rand(ry(i), 1)*(ry(i-1)*n(i-1))));
                rkick = numel(ind2);
                curbl = zeros(rkick*n(i)*ry(i+1), nx);
                for j=1:nx
                    cr1 = reshape(crX{i,j}, rx(i,j)*n(i), rx(i+1,j));
                    cr1 = cr1*Rx{i+1,j};
                    cr1 = reshape(cr1, rx(i,j), n(i)*ry(i+1));
                    cr2 = reshape(crX{i-1,j}, rx(i-1,j), n(i-1)*rx(i,j));
                    cr2 = Rx{i-1,j}*cr2; % now its size rx
                    cr2 = reshape(cr2, ry(i-1)*n(i-1), rx(i,j));
                    cr2 = cr2(ind2, :);
                    curbl(:,j) = reshape(cr2*cr1, rkick*n(i)*ry(i+1), 1);
                end;
                % Call the function
                uk = funs(curbl); % rnr x d2
                uk = reshape(uk, rkick*n(i)*ry(i+1), d2);
                uk = reshape(uk.', d2*rkick*n(i), ry(i+1));
                uk = uk / Ry{i+1};
                uk = reshape(uk, d2*rkick, n(i)*ry(i+1));
                if (strcmp(pcatype, 'svd'))
                    [vk,sk,uk]=svd(uk, 'econ');
                    uk = uk(:,1:min(kickrank, size(uk,2)));
                else
                    uk = uchol(uk, kickrank+1);
                    uk = uk(:,end:-1:max(end-kickrank+1,1));
                end;
            else
                uk = rand(n(i)*ry(i+1), kickrank);
            end;
%             uk = rand(n(i)*ry(i+1), kickrank);
            [v,rv]=qr([v,uk], 0);
            radd = size(uk,2);
        end;
        u = [u, zeros(d2*ry(i), radd)];
        u = u*(rv.');
        r = size(v,2);

        cr2 = cry{i-1};
        cr2 = reshape(cr2, ry(i-1)*n(i-1), ry(i));
        u = reshape(u, d2, ry(i)*r);
        u = reshape(u.', ry(i), r*d2);
        u = cr2*u;

        u = reshape(u, ry(i-1)*n(i-1)*r, d2);
        u = reshape(u.', d2, ry(i-1), n(i-1), r);
        v = reshape(v.', r, n(i), ry(i+1));

        % Stuff back
        ry(i) = r;
        cry{i-1} = u;
        cry{i} = v;

        % Recompute left interface matrices
        % Interface matrix for Y
        Ry{i} = reshape(v, ry(i)*n(i), ry(i+1))*Ry{i+1};
        Ry{i} = reshape(Ry{i}, ry(i), n(i)*ry(i+1));
        curind = maxvol2(Ry{i}.','qr',do_qr);
        Ry{i} = Ry{i}(:, curind);
        % Interface matrices for X
        for j=1:nx
            Rx{i,j} = reshape(crX{i,j}, rx(i,j)*n(i), rx(i+1,j));
            Rx{i,j} = Rx{i,j}*Rx{i+1,j};
            Rx{i,j} = reshape(Rx{i,j}, rx(i,j), n(i)*ry(i+1));
            Rx{i,j} = Rx{i,j}(:, curind);
        end;
    elseif ((dir>0)&&(i==d))
        % Just stuff back the last core
        newy = u(:,1:r)*diag(s(1:r))*(v(:,1:r)');
        newy = reshape(newy, ry(i)*n(i)*ry(i+1), d2);
        cry{i} = reshape(newy.', d2, ry(i), n(i), ry(i+1));
    elseif ((dir<0)&&(i==1))
        % Just stuff back the last core
        newy = u(:,1:r)*diag(s(1:r))*(v(:,1:r)');
        newy = reshape(newy, d2, ry(i), n(i), ry(i+1));
        cry{i} = newy;
    end;

    i = i+dir;
    % Reversing, residue check, etc
    cur_order(order_index) = cur_order(order_index) - dir;
    % New direction
    if (cur_order(order_index)==0)
        order_index = order_index+1;

        if (verb>0)
            fprintf('=multifuncrs= sweep %d{%d}, max_dy: %3.3e, erank: %g\n', swp, order_index-1, max_dy, sqrt(ry(1:d)'*(n.*ry(2:d+1))/sum(n)));
        end;

        if (last_sweep)
            break;
        end;

        if (max_dy<eps)&&(dir<0)
            last_sweep=true;
            kickrank=0;
        end;

        if (order_index>numel(cur_order)) % New global sweep
            cur_order = block_order;
            order_index = 1;
            %residue
            if (last_sweep)
                cur_order = d-1;
            end;

            max_dy = 0;
            swp = swp+1;
        end;

        dir = sign(cur_order(order_index));
        i = i+dir;
    end;
end

cry{d} = permute(cry{d}, [2,3,1]); % d2 is r(d+1)
y = cell2core(y, cry);

end

author: oseledets

repo: TT-Toolbox


text:
function outcome=getOutcomeStr(obj, verbosity)
% give string indicating the outcome of the test
%
% outcome=getOutcomeStr(obj, verbosity)
%
% Inputs:
%   obj                 MOxUnitTestOutcome object
%   verbosity           Integer in the range [0,1,2]
%
% Returns:
%   outcome:            Depending on the value of verbosity
%                       0: ''
%                       1: 's'
%                       2: 'skipped'
%

    outcome_cell={'','s','skipped'};
    outcome=outcome_cell{verbosity+1};

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function dq = rotvel2dquat(omega,axis)

% ROTVEL2DQUAT  constructs a rotational velocity dual quaternion
%
%    DQ = ROTVEL2DQUAT(OMEGA,AXIS) returns the dual quaternion DQ which is
%     the dual quaternion representation of the rotational velocity of axis
%     AXIS and of angular speed OMEGA [deg/s].
%       -   OMEGA represents the angular speed [deg/s]. It is a scalar or a
%               N-vector (the i-th element represents the angular speed of
%               the i-th rotational velocity) where N is the number of
%               rotational velocities.
%       -   AXIS represents the rotational velocity axis. It is a 3-vector
%               or a 3*N array (column i represents the rotation axis of
%               the i-th rotational velocity).
%       -   DQ is a 8*N array representing the rotational velocity dual
%              quaternion [rad/s]
%               DQ = [0 axis'*omega 0 0 0 0]
%
% See also ROT2DQUAT, TRANSVEL2DQUAT

somega = size(omega);
saxis = size(axis);
if saxis == [1 3], axis = axis'; saxis = size(axis); end
if somega(1) > 1, omega = omega'; somega = size(omega); end

% wrong size
if somega(1) ~= 1 || saxis(1) ~= 3
    error('DualQuaternion:rotvel2dquat:wrongsize',...
        '%d rows in array omega and %d rows in array axis. It should be 1 for omega and 3 for axis.',...
        somega(1),saxis(1));
end

nomega = somega(2);
naxis = saxis(2);
if nomega ~= naxis
    error('DualQuaternion:rotvel2dquat:wrongFormat',...
        '%d elements in array omega and %d axis in array axis. It should be the same number for both.',...
        nomega,naxis);
end

% normalization of the axis vector (if necessary)
n = naxis;
n2 = sum(axis.^2).^0.5;

indnot0 = find(n2>0);
axis(:,indnot0) = axis(:,indnot0)./repmat(n2(:,indnot0),3,1);

% construction of the rotation dual quaternion
dq = zeros(8,n);
omegatab = repmat(omega,3,1)*pi/180;
dq(2:4,:) = omegatab.*axis;


text:
function is=lpcrf2is(rf)
%LPCRF2IS Convert reflection coefficients to inverse sines IS=(RF)

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpcrf2is.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

is=asin(rf)*2/pi;

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
%CODINGCONVENTIONS specify coding conventions for MALIMPA m-files
%
%   General:
%
%   This file specifies coding convention for future files. Some files may
%   not necessarily follows the requirements, as they were created before
%   conventions.
%
%   Name of functions:
%
%   Java convention for file names must be used:
%   - Each function name must be written in lower case, using english
%       denomination.
%   - For compound function names, the first letter of each component after
%       the first one is capitalized.
%       Ex: 'distancePointPlane'.
%
%   Other convention should be used to avoid name conflicts:
%   - For function used as image filters, prefix 'im' may be used. Ex :
%       'imDirFilter'.
%   - For functions which work on a complex data type, the name of the data
%       type should be used as a prefix of the function name.
%       Ex: functions working on polygons have a prototype on the form
%       'polygonLength', 'polygonArea'...
%   - Functions working with geometrical primitives can use a suffix for
%       indicating the dimension of application.
%       Ex: 'distancePointLine3d'.
%       This convention was not used for the beginning of the library
%       development, but new functions will try to conform to it.
%
%   Function header:
%
%   - First line of documentation must contain the name of the function,
%   followed by a short description. Name and function may be separated by
%   colon. Name must be capitalized, without punctuation before
%   description (to remain valid with matlab documenbtation tools).
%   Ex:   'SLICER Interactive visualization of 3D images'
%
%   - Main part of header describes usage of function/script with various
%   parameters. Function name should be the real name. Variables names
%   should be capitalized. This does not follow matlab's standard, but
%   allows to copy-paste pieces of code.
%   Ex: 'RES = myBasicFunction(PARAM1, PARAM2, OPTION);'
%
%   - The description of each calling syntax should be written using
%   conjugated verbs.
%
%   - End of header should contain name of author, with date of creation.
%   It may also contains author affiliation (company ...), and date of last
%   revision.
%
%
%   Library Organization:
%
%   The principle is to gather functions which share a similar domain of
%   application. Ex: image processing filters, geometry functions...
%   It seems better to have one level of complexity for the directory
%   structure.
%
%   However, it could be clearer to create specific subdirectories for
%   functions which operate on complex data (such as polygons, or special
%   types of curves). This allows to add a Contents file which describes
%   the model used by functions.
%
%
%   Function code:
%
%   - Please comment the code ! For algorithmic parts, as much comments as
%   code is a good start ...
%   - Indentation level should be marked by four space characters.
%
%
%   Tips & Tricks:
%
%   - The use of '...' at the end of a line can cause compatibility
%   problems from Linux to Windows.
%   - matlab help uses the file 'Contents.m', not 'contents.m' !
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - MIA MathCell
%   created the 23/03/2005
%

%   HISTORY
%   05/03/2007 add convention for datatype prefix for function names, and
%       library organization

help codingConventions;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

function prep_deploy()
    % to be run to set the directory up for deployment
    u = userpath;
    matlab_path = u(1:end-1);

    code_root_dir = fullfile(fileparts(which(mfilename)), '..');
    extern_src_rel_dst = 'extern_src';
    algos_dir_dst = fullfile(code_root_dir, extern_src_rel_dst);
    utils_dir_dst = fullfile(code_root_dir, 'utils');
    algos_dir_src = fullfile(code_root_dir, '..', 'extern_src');

    % delete old directories & files (as a result of previous run of this script)
    if_exist_del(fullfile(code_root_dir, 'data'));
    if_exist_del(fullfile(algos_dir_dst, 'fuxin_lib_src'));
    if_exist_del(fullfile(algos_dir_dst, 'segmentation'));
    if_exist_del(fullfile(algos_dir_dst, 'stein_boundaryprocessing'));
    if_exist_del(fullfile(algos_dir_dst, 'toolboxes'));
    if_exist_del(fullfile(algos_dir_dst, 'utils'));
    if_exist_del(fullfile(utils_dir_dst, 'drawFigFrames.m'));
    if_exist_del(fullfile(utils_dir_dst, 'timerTicks.m'));
    if_exist_del(fullfile(utils_dir_dst, 'compute_error_metric.m'));
    if_exist_del(fullfile(utils_dir_dst, 'duplicateElems.m'));

    % copy files needed from matlab user directory
    copyfile(fullfile(matlab_path, 'drawFigFrames.m'), utils_dir_dst);
    copyfile(fullfile(matlab_path, 'timerTicks.m'), utils_dir_dst);
    copyfile(fullfile(matlab_path, 'duplicateElems.m'), utils_dir_dst);

    % copy necessary files from fuxin library
    mkdir(fullfile(algos_dir_dst, 'fuxin_lib_src'));
    copyfile(fullfile(code_root_dir, '..', 'fuxin_lib_src', 'boosting'), fullfile(algos_dir_dst, 'fuxin_lib_src', 'boosting'));
    copyfile(fullfile(code_root_dir, '..', 'fuxin_lib_src', 'myqueue_1.1'), fullfile(algos_dir_dst, 'fuxin_lib_src', 'myqueue_1.1'));
    copyfile(fullfile(code_root_dir, '..', 'fuxin_lib_src', '@classregtree_fuxin'), fullfile(algos_dir_dst, 'fuxin_lib_src', '@classregtree_fuxin'));

    % copy stein boundary processing
    copyfile(fullfile(algos_dir_src, 'segmentation', 'stein_boundaryprocessing'), fullfile(algos_dir_dst, 'stein_boundaryprocessing'));

    copyfile(fullfile(algos_dir_src, 'compute_error_metric.m'), utils_dir_dst);
end

function if_exist_del(delpath)
    if exist(delpath, 'dir')
        rmdir(delpath, 's');
    elseif exist(delpath, 'file')
        delete(delpath);
    end
end

author: Cloud-CV

repo: object-proposals


text:
function [ output ] = tapas_rdcm_compute_signals(DCM, output, options)
% [ output ] = tapas_rdcm_compute_signals(DCM, output, options)
%
% Computes true and predicted signals
%
%   Input:
%   	DCM         - model structure
%       output      - model inversion results
%       options     - estimation options
%
%   Output:
%   	output      - model inversion results with model fits
%

% ----------------------------------------------------------------------
%
% Authors: Stefan Fraessle (stefanf@biomed.ee.ethz.ch), Ekaterina I. Lomakina
%
% Copyright (C) 2016-2022 Translational Neuromodeling Unit
%                         Institute for Biomedical Engineering
%                         University of Zurich & ETH Zurich
%
% This file is part of the TAPAS rDCM Toolbox, which is released under the
% terms of the GNU General Public License (GPL), version 3.0 or later. You
% can redistribute and/or modify the code under the terms of the GPL. For
% further see COPYING or <http://www.gnu.org/licenses/>.
%
% Please note that this toolbox is in an early stage of development. Changes
% are likely to occur in future releases.
%
% ----------------------------------------------------------------------

% true or measured signal
output.signal.y_source = DCM.Y.y(:);

% true (deterministic) signal / VBL signal
if strcmp(options.type,'s')

    % if signal should be computed
    if ( options.compute_signal )

        % noise-free signal
        DCMs = tapas_rdcm_generate(DCM, options, Inf);
        output.signal.y_clean = DCMs.Y.y(:);

        % compute the MSE of the noisy data
        output.residuals.y_mse_clean = mean((output.signal.y_source - output.signal.y_clean).^2);

    end

else

    % VBL predicted signal
    if ( isfield(DCM,'y') )
        output.signal.y_pred_vl     = DCM.y(:);
        output.residuals.y_mse_vl   = mean((output.signal.y_source - output.signal.y_pred_vl).^2);
    else
        output.residuals.y_mse_vl   = [];
    end
end

% store true or measured temporal derivative (in frequency domain)
yd_source_fft                           = output.temp.yd_source_fft;
yd_source_fft(~isfinite(yd_source_fft)) = 0;
output.signal.yd_source_fft             = yd_source_fft(:);

% if signal in time domain should be computed
if ( options.compute_signal )

    % get rDCM parameters
    DCM.Tp   = output.Ep;

    % no baseline for simulated data
    if ( strcmp(options.type,'s') )
        DCM.Tp.baseline = zeros(size(DCM.Tp.C,1),1);
    end

    % increase sampling rate
    r_dt     = DCM.Y.dt/DCM.U.dt;
    DCM.Y.dt = DCM.U.dt;

    % posterior probability (for sparse rDCM)
    if ( isfield(output,'Ip') )
        DCM.Tp.A = DCM.Tp.A .* output.Ip.A;
        DCM.Tp.C = DCM.Tp.C .* output.Ip.C;
    end

    % generate predicted signal (tapas_rdcm_generate)
    DCM_rDCM = tapas_rdcm_generate(DCM, options, Inf);

    % add the confounds to predicted time series
    for t = 1:size(DCM.Y.y,1)
        for r = 1:size(DCM.Y.y,2)
            DCM_rDCM.Y.y(t,r) = DCM_rDCM.Y.y(t,r) + DCM_rDCM.Tp.baseline(r,:) * DCM_rDCM.U.X0((1+r_dt*(t-1)),:)';
        end
    end

    % turn into vector
    output.signal.y_pred_rdcm = DCM_rDCM.Y.y(:);

end

% store predicted temporal derivative (in frequency domain)
yd_pred_rdcm_fft                = output.temp.yd_pred_rdcm_fft;
output.signal.yd_pred_rdcm_fft  = yd_pred_rdcm_fft(:);

% remove the temp structure
output = rmfield(output,'temp');

% asign the region names
if ( isfield(DCM.Y,'name') )
    output.signal.name = DCM.Y.name;
end

% compute the MSE of predicted signal
if ( options.compute_signal )
    output.residuals.y_mse_rdcm     = mean((output.signal.y_source - output.signal.y_pred_rdcm).^2);
    output.residuals.R_rdcm         = output.signal.y_source - output.signal.y_pred_rdcm;
end

% store the driving inputs
output.inputs.u     = DCM.U.u;
output.inputs.name  = DCM.U.name;
output.inputs.dt    = DCM.U.dt;

end

author: translationalneuromodeling

repo: tapas


text:
classdef TotalCoordinatesCalculator < handle

    properties (Access = public, Abstract)
        c
        theta
        nodes
        vertCoord
        boundCoord
        div
        totalCoord
    end

    methods (Access = public, Static)

        function obj = create(cParams)
            obj = TotalCoordinatesCalculatorFactory.create(cParams);
        end

    end

    methods (Access = protected)

        function init(obj,cParams)
            obj.c  = cParams.c;
            obj.theta = cParams.theta;
            obj.nodes = cParams.nodes;
            obj.vertCoord = cParams.vertCoord;
            obj.boundCoord = cParams.boundCoord;
            obj.div = cParams.div;
            obj.totalCoord = [];
        end

        function initBoundary(obj)
            boundNodes = obj.nodes.bound;
            obj.totalCoord(1:boundNodes,:) = obj.boundCoord;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% ==============================================================================
%
%    Software License Agreement (BSD License)
%    Copyright (c) 2019
%    (www.aimlab.wpi.edu)
%
%    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 authors 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.
%
%    \author:    <http://www.aimlab.wpi.edu>
%    \author:    <amunawar@wpi.edu>
%    \author:    Adnan Munawar
%    \version:   0.1$
% ==============================================================================

function [its,sizePath,run_time] =  RRTextend3D(dim,segmentLength,random_world,show_output)
% dim = 2;
% radius =0;
% segmentLength = 5;
% random_world = 0;
% n_its = 1000;
% standard length of path segments
if dim ==2
    start_cord = [5,5];
    goal_cord = [95,95];

else

    start_cord = [5,5,5];
    goal_cord = [95,95,95];
end

% create random world
Size = 100;
NumObstacles = 100;

if random_world ==1
    world = createWorld(NumObstacles,ones(1,dim)*Size,zeros(1,dim),dim);
else
    [world NumObstacles] = createKnownWorld(ones(1,dim)*Size,[0;0;0],dim);
end
% randomly select start and end nodes
%start_node = generateRandomNode(world,dim)
%end_node   = generateRandomNode(world,dim)
start_node = [start_cord,0,0,0];
end_node = [goal_cord,0,0,0];
% establish tree starting with the start node
tree = start_node;

a = clock;

% check to see if start_node connects directly to end_node
if ( (norm(start_node(1:dim)-end_node(1:dim))<segmentLength )...
        &&(collision(start_node,end_node,world,dim)==0) )
    path = [start_node; end_node];
else

    its = 0;
    numPaths = 0;
    flag = 0;
    while numPaths < 1,
        [tree,flag] = extendTree(tree,end_node,segmentLength,world,flag,dim);
        numPaths = numPaths + flag;
        its = its+1;
    end

end

% find path with minimum cost to end_node
path = findMinimumPath(tree,end_node,dim);
sizePath = size(path,1);

b = clock;
run_time = 3600*(b(4)-a(4)) + 60 * (b(5)-a(5)) + (b(6) - a(6));

if show_output == 1
    figure;
    plotExpandedTree(world,tree,dim);
    plotWorld(world,path,dim);
    % figure(2);
    % plotWorld(world,path,dim);
    %plotExpandedTree(world,tree,dim);
end
end

function world = createWorld(NumObstacles, endcorner, origincorner,dim)

if dim == 2

    % check to make sure that the region is nonempty
    if (endcorner(1) <= origincorner(1)) | (endcorner(2) <= origincorner(2))
        disp('Not valid corner specifications!')
        world=[];

        % create world data structure
    else
        world.NumObstacles = NumObstacles;
        world.endcorner = endcorner;
        world.origincorner = origincorner;

        % create NumObstacles
        maxRadius = min(endcorner(1)- origincorner(1), endcorner(2)-origincorner(2));
        maxRadius = 5*maxRadius/NumObstacles/2;
        for i=1:NumObstacles,
            % randomly pick radius
            world.radius(i) = maxRadius*rand;
            % randomly pick center of obstacles
            cx = origincorner(1) + world.radius(i)...
                + (endcorner(1)-origincorner(1)-2*world.radius(i))*rand;
            cy = origincorner(2) + world.radius(i)...
                + (endcorner(2)-origincorner(2)-2*world.radius(i))*rand;
            world.cx(i) = cx;
            world.cy(i) = cy;
        end
    end

elseif dim ==3;
    % check to make sure that the region is nonempty
    if (endcorner(1) <= origincorner(1)) || (endcorner(2) <= origincorner(2)) || (endcorner(3) <= origincorner(3))
        disp('Not valid corner specifications!')
        world=[];

        % create world data structure
    else
        world.NumObstacles = NumObstacles;
        world.endcorner = endcorner;
        world.origincorner = origincorner;

        % create NumObstacles
        bounds = [endcorner(1)- origincorner(1), endcorner(2)-origincorner(2), endcorner(3)-origincorner(3)];
        maxRadius = min(bounds);
        maxRadius = 5*maxRadius/NumObstacles;
        for i=1:NumObstacles,
            % randomly pick radius
            world.radius(i) = maxRadius*rand;
            % randomly pick center of obstacles
            cx = origincorner(1) + world.radius(i)...
                + (endcorner(1)-origincorner(1)-2*world.radius(i))*rand;
            cy = origincorner(2) + world.radius(i)...
                + (endcorner(2)-origincorner(2)-2*world.radius(i))*rand;
            cz = origincorner(2) + world.radius(i)...
                + (endcorner(2)-origincorner(2)-2*world.radius(i))*rand;
            world.cx(i) = cx;
            world.cy(i) = cy;
            world.cz(i) = cz;
        end
    end
end
end

function [world NumObstacles] = createKnownWorld(endcorner, origincorner,dim)
NumObstacles = 5;
if dim == 2
    % check to make sure that the region is nonempty
    if (endcorner(1) <= origincorner(1)) | (endcorner(2) <= origincorner(2)),
        disp('Not valid corner specifications!')
        world=[];
        % create world data structure
    else
        world.NumObstacles = NumObstacles;
        world.endcorner = endcorner;
        world.origincorner = origincorner;

        % create NumObstacles
        maxRadius = 10;

        world.radius(1) = maxRadius;
        cx = 50;
        cy = 50;
        world.cx(1) = cx;
        world.cy(1) = cy;

        world.radius(2) = maxRadius;
        cx = 75;
        cy = 25;
        world.cx(2) = cx;
        world.cy(2) = cy;

        world.radius(3) = maxRadius;
        cx = 25;
        cy = 75;
        world.cx(3) = cx;
        world.cy(3) = cy;

        world.radius(4) = maxRadius;
        cx = 25;
        cy = 25;
        world.cx(4) = cx;
        world.cy(4) = cy;

        world.radius(5) = maxRadius;
        cx = 75;
        cy = 75;
        world.cx(5) = cx;
        world.cy(5) = cy;
    end

elseif dim == 3

    NumObstacles = 9;
    % check to make sure that the region is nonempty
    if (endcorner(1) <= origincorner(1)) | (endcorner(2) <= origincorner(2)) | (endcorner(3) <= origincorner(3)),
        disp('Not valid corner specifications!')
        world=[];

        % create world data structure
    else
        world.NumObstacles = NumObstacles;
        world.endcorner = endcorner;
        world.origincorner = origincorner;

        % create NumObstacles
        maxRadius = 10;

        world.radius(1) = maxRadius;
        cx = 50;
        cy = 50;
        cz = 50;
        world.cx(1) = cx;
        world.cy(1) = cy;
        world.cz(1) = cz;

        world.radius(2) = maxRadius;
        cx = 25;
        cy = 25;
        cz = 25;
        world.cx(2) = cx;
        world.cy(2) = cy;
        world.cz(2) = cz;

        world.radius(3) = maxRadius;
        cx = 75;
        cy = 75;
        cz = 75;
        world.cx(3) = cx;
        world.cy(3) = cy;
        world.cz(3) = cz;

        world.radius(4) = maxRadius;
        cx = 25;
        cy = 25;
        cz = 75;
        world.cx(4) = cx;
        world.cy(4) = cy;
        world.cz(4) = cz;

        world.radius(5) = maxRadius;
        cx = 75;
        cy = 75;
        cz = 25;
        world.cx(5) = cx;
        world.cy(5) = cy;
        world.cz(5) = cz;

        world.radius(6) = maxRadius;
        cx = 25;
        cy = 75;
        cz = 25;
        world.cx(6) = cx;
        world.cy(6) = cy;
        world.cz(6) = cz;

        world.radius(7) = maxRadius;
        cx = 75;
        cy = 25;
        cz = 25;
        world.cx(7) = cx;
        world.cy(7) = cy;
        world.cz(7) = cz;

        world.radius(8) = maxRadius;
        cx = 75;
        cy = 25;
        cz = 75;
        world.cx(8) = cx;
        world.cy(8) = cy;
        world.cz(8) = cz;

        world.radius(9) = maxRadius;
        cx = 25;
        cy = 75;
        cz = 75;
        world.cx(9) = cx;
        world.cy(9) = cy;
        world.cz(9) = cz;
    end
end
end

function node=generateRandomNode(world,dim)

if dim ==2;
    % randomly pick configuration
    px       = (world.endcorner(1)-world.origincorner(1))*rand;
    py       = (world.endcorner(2)-world.origincorner(2))*rand;

    chi      = 0;
    cost     = 0;
    node     = [px, py, chi, cost, 0];

    % check collision with obstacle
    while collision(node, node, world,dim),
        px       = (world.endcorner(1)-world.origincorner(1))*rand;
        py       = (world.endcorner(2)-world.origincorner(2))*rand;

        chi      = 0;
        cost     = 0;
        node     = [px, py, chi, cost, 0];
    end

elseif dim ==3;
    % randomly pick configuration
    px       = (world.endcorner(1)-world.origincorner(1))*rand;
    py       = (world.endcorner(2)-world.origincorner(2))*rand;
    pz       = (world.endcorner(3)-world.origincorner(3))*rand;

    chi      = 0;
    cost     = 0;
    node     = [px, py, pz, chi, cost, 0];

    % check collision with obstacle
    while collision(node, node, world,dim),
        px       = (world.endcorner(1)-world.origincorner(1))*rand;
        py       = (world.endcorner(2)-world.origincorner(2))*rand;
        pz       = (world.endcorner(3)-world.origincorner(3))*rand;

        chi      = 0;
        cost     = 0;
        node     = [px, py, pz, chi, cost, 0];
    end

end

end

function collision_flag = collision(node, parent, world,dim)

collision_flag = 0;

for i=1:dim
    if (node(i)>world.endcorner(i))|(node(i)<world.origincorner(i))
        collision_flag = 1;
    end
end

if collision_flag == 0 && dim ==2
    for sigma = 0:.2:1,
        p = sigma*node(1:dim) + (1-sigma)*parent(1:dim);
        % check each obstacle
        for i=1:world.NumObstacles,
            if (norm([p(1);p(2)]-[world.cx(i); world.cy(i)])<=1*world.radius(i)),
                collision_flag = 1;
                break;
            end
        end
    end

elseif collision_flag == 0 && dim ==3
    for sigma = 0:.2:1,
        p = sigma*node(1:dim) + (1-sigma)*parent(1:dim);
        % check each obstacle
        for i=1:world.NumObstacles,
            if (norm([p(1);p(2);p(3)]-[world.cx(i); world.cy(i); world.cz(i)])<=1*world.radius(i)),
                collision_flag = 1;
                break;
            end
        end
    end
end
end

function flag = canEndConnectToTree(tree,end_node,minDist,world,dim)
flag = 0;
% check only last node added to tree since others have been checked
if ( (norm(tree(end,1:dim)-end_node(1:dim))<minDist)...
        & (collision(tree(end,1:dim), end_node(1:dim), world,dim)==0) ),
    flag = 1;
end

end

function [new_tree,flag] = extendTree(tree,end_node,segmentLength,world,flag_chk,dim)
r = 0;
flag1 = 0;
while flag1==0,
    % select a random point
    randomPoint = ones(1,dim);
    for i=1:dim
        randomPoint(1,i) = (world.endcorner(i)-world.origincorner(i))*rand;
    end

    % find leaf on node that is closest to randomPoint
    tmp = tree(:,1:dim)-ones(size(tree,1),1)*randomPoint;
    sqrd_dist = sqr_eucl_dist(tmp,dim);
    [min_dist,idx] = min(sqrd_dist);
    min_parent_idx = idx;

    new_point = (randomPoint-tree(idx,1:dim));
    new_point = tree(idx,1:dim)+(new_point/norm(new_point))*segmentLength;

    min_cost  = cost_np(tree(idx,:),new_point,dim);
    new_node  = [new_point, 0, min_cost, idx];

    if collision(new_node, tree(idx,:), world,dim)==0

        tmp_dist = tree(:,1:dim)-(ones(size(tree,1),1)*new_point);
        dist = sqr_eucl_dist(tmp_dist,dim);
        near_idx = find(dist <= r^2);

        if size(near_idx,1)>1
            size_near = size(near_idx,1);

            for i = 1:size_near
                if collision(new_node, tree(near_idx(i),:), world,dim)==0

                    cost_near = tree(near_idx(i),dim+2)+line_cost(tree(near_idx(i),:),new_point,dim);

                    if  cost_near < min_cost
                        min_cost = cost_near;
                        min_parent_idx = near_idx(i);
                    end

                end
            end
        end

        new_node = [new_point, 0 , min_cost, min_parent_idx];
        new_tree = [tree; new_node];
        new_node_idx = size(new_tree,1);

        if size(near_idx,1)>1
            reduced_idx = near_idx;
            for j = 1:size(reduced_idx,1)
                near_cost = new_tree(reduced_idx(j),dim+2);
                lcost = line_cost(new_tree(reduced_idx(j),:),new_point,dim);
                if near_cost > min_cost + lcost ...
                        && collision(new_tree(reduced_idx(j),:),new_node,world,dim)
                    before = new_tree(reduced_idx(j),dim+3)
                    new_tree(reduced_idx(j),dim+3) = new_node_idx;
                    after = new_tree(reduced_idx(j),dim+3)
                end

            end
        end
        flag1=1;
    end
end

if flag_chk == 0
    % check to see if new node connects directly to end_node
    if ( (norm(new_node(1:dim)-end_node(1:dim))<segmentLength )...
            && (collision(new_node,end_node,world,dim)==0) )
        flag = 1;
        new_tree(end,dim+1)=1;  % mark node as connecting to end.
    else
        flag = 0;
    end

else flag = 1;
end
end

function e_dist = sqr_eucl_dist(array,dim)

sqr_e_dist = zeros(size(array,1),dim);
for i=1:dim

    sqr_e_dist(:,i) = array(:,i).*array(:,i);

end
e_dist = zeros(size(array,1),1);
for i=1:dim

    e_dist = e_dist+sqr_e_dist(:,i);

end

end

%calculate the cost from a node to a point
function [cost] = cost_np(from_node,to_point,dim)

diff = from_node(:,1:dim) - to_point;
eucl_dist = norm(diff);
cost = from_node(:,dim+2) + eucl_dist;

end

%calculate the cost from a node to a node
function [cost] = cost_nn(from_node,to_node,dim)

diff = from_node(:,1:dim) - to_node(:,1:dim);
eucl_dist = norm(diff);
cost = from_node(:,dim+2) + eucl_dist;

end

function [cost] = line_cost(from_node,to_point,dim)
diff = from_node(:,1:dim) - to_point;
cost = norm(diff);
end

function path = findMinimumPath(tree,end_node,dim)

% find nodes that connect to end_node
connectingNodes = [];
for i=1:size(tree,1),
    if tree(i,dim+1)==1,
        connectingNodes = [connectingNodes ; tree(i,:)];
    end
end

% find minimum cost last node
[tmp,idx] = min(connectingNodes(:,dim+2));

% construct lowest cost path
path = [connectingNodes(idx,:); end_node];
parent_node = connectingNodes(idx,dim+3);
while parent_node>1,
    parent_node = tree(parent_node,dim+3);
    path = [tree(parent_node,:); path];
end

end

function plotExpandedTree(world,tree,dim)
ind = size(tree,1);
while ind>0
    branch = [];
    node = tree(ind,:);
    branch = [ branch ; node ];
    parent_node = node(dim+3);
    while parent_node > 1
        cur_parent = parent_node;
        branch = [branch; tree(parent_node,:)];
        parent_node = tree(parent_node,dim+3);
    end
    ind = ind - 1;

    if dim == 2
        X = branch(:,1);
        Y = branch(:,2);

        p = plot(X,Y);
        set(p,'Color','r','LineWidth',0.5,'Marker','.','MarkerEdgeColor','g');
        hold on;

    elseif dim == 3
        X = branch(:,1);
        Y = branch(:,2);
        Z = branch(:,3);

        p = plot3(X,Y,Z);
        set(p,'Color','r','LineWidth',0.5,'Marker','.','MarkerEdgeColor','g');
        hold on;
    end
end
end

function plotWorld(world,path,dim)
% the first element is the north coordinate
% the second element is the south coordinate
if dim ==2

    N = 10;
    th = 0:2*pi/N:2*pi;
    axis([world.origincorner(1),world.endcorner(1),...
        world.origincorner(2), world.endcorner(2)]);
    hold on

    for i=1:world.NumObstacles,
        X = world.radius(i)*sin(th) + world.cx(i);
        Y = world.radius(i)*cos(th) + world.cy(i);
        fill(X,Y,'blue');
    end

    X = path(:,1);
    Y = path(:,2);
    p = plot(X,Y);

elseif dim ==3
    axis([world.origincorner(1),world.endcorner(1),...
        world.origincorner(2), world.endcorner(2),...
        world.origincorner(3), world.endcorner(3)]);
    hold on

    for i=1:world.NumObstacles,
        [X Y Z] = sphere(10);
        X = (X*world.radius(i));
        Y = (Y*world.radius(i));
        Z = (Z*world.radius(i));
        surf(X+world.cx(i),Y+world.cy(i),Z+world.cz(i));
        colormap([0.5 0.2 0.3]);
    end

    X = path(:,1);
    Y = path(:,2);
    Z = path(:,3);
    p = plot3(X,Y,Z);
end
set(p,'Color','black','LineWidth',3)
xlabel('X axis');
ylabel('Y axis');
zlabel('Z axis');
title('RRT Extend Algorithm');
end

author: adnanmunawar

repo: matlab-rrt-variants


text:
% eegplugin_dipfit() - DIPFIT plugin version 2.0 for EEGLAB menu.
%                      DIPFIT is the dipole fitting Matlab Toolbox of
%                      Robert Oostenveld (in collaboration with A. Delorme).
%
% Usage:
%   >> eegplugin_dipfit(fig, trystrs, catchstrs);
%
% Inputs:
%   fig        - [integer] eeglab figure.
%   trystrs    - [struct] "try" strings for menu callbacks.
%   catchstrs  - [struct] "catch" strings for menu callbacks.
%
% Notes:
%   To create a new plugin, simply create a file beginning with "eegplugin_"
%   and place it in your eeglab folder. It will then be automatically
%   detected by eeglab. See also this source code internal comments.
%   For eeglab to return errors and add the function's results to
%   the eeglab history, menu callback must be nested into "try" and
%   a "catch" strings. For more information on how to create eeglab
%   plugins, see http://www.sccn.ucsd.edu/eeglab/contrib.html
%
% Author: Arnaud Delorme, CNL / Salk Institute, 22 February 2003
%
% See also: eeglab()

% Copyright (C) 2003 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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-1.07  USA

function vers = eegplugin_dipfit(fig, trystrs, catchstrs)

    vers = 'dipfit2.2';
    if nargin < 3
        error('eegplugin_dipfit requires 3 arguments');
    end;

    % add dipfit folder to path
    % -----------------------
    if ~exist('ft_dipolefitting')
        p = which('eegplugin_dipfit');
        disp('Fieldtrip functions for dipole localization not found, removing Dipfit2');
        p = p(1:findstr(p,'eegplugin_dipfit.m')-1);
        rmpath(p);
    end;

    % find tools menu
    % ---------------
    menu = findobj(fig, 'tag', 'tools');
    % tag can be
    % 'import data'  -> File > import data menu
    % 'import epoch' -> File > import epoch menu
    % 'import event' -> File > import event menu
    % 'export'       -> File > export
    % 'tools'        -> tools menu
    % 'plot'         -> plot menu

    % command to check that the '.source' is present in the EEG structure
    % -------------------------------------------------------------------
    check_dipfit = [trystrs.no_check 'if ~isfield(EEG, ''dipfit''), error(''Run the dipole setting first''); end;'  ...
                    'if isempty(EEG.dipfit), error(''Run the dipole setting first''); end;'  ];
    check_dipfitnocheck = [ trystrs.no_check 'if ~isfield(EEG, ''dipfit''), error(''Run the dipole setting first''); end; ' ];
    check_chans = [ '[EEG tmpres] = eeg_checkset(EEG, ''chanlocs_homogeneous'');' ...
                       'if ~isempty(tmpres), h(tmpres), end; clear tmpres;' ];

    % menu callback commands
    % ----------------------
    comsetting = [ trystrs.check_ica check_chans '[EEG LASTCOM]=pop_dipfit_settings(EEG);'    catchstrs.store_and_hist ];
    combatch   = [ check_dipfit check_chans  '[EEG LASTCOM] = pop_dipfit_gridsearch(EEG);'    catchstrs.store_and_hist ];
    comfit     = [ check_dipfitnocheck check_chans [ 'EEG = pop_dipfit_nonlinear(EEG); ' ...
                        'LASTCOM = ''% === History not supported for manual dipole fitting ==='';' ]  catchstrs.store_and_hist ];
    comauto    = [ check_dipfit check_chans  '[EEG LASTCOM] = pop_multifit(EEG);'        catchstrs.store_and_hist ];
    % preserve the '=" sign in the comment above: it is used by EEGLAB to detect appropriate LASTCOM
    complot    = [ check_dipfit check_chans 'LASTCOM = pop_dipplot(EEG);'                     catchstrs.add_to_hist ];

    % create menus
    % ------------
    submenu = uimenu( menu, 'Label', 'Locate dipoles using DIPFIT 2.x', 'separator', 'on');
    uimenu( submenu, 'Label', 'Head model and settings'  , 'CallBack', comsetting);
    uimenu( submenu, 'Label', 'Coarse fit (grid scan)'   , 'CallBack', combatch);
    uimenu( submenu, 'Label', 'Fine fit (iterative)'     , 'CallBack', comfit);
    uimenu( submenu, 'Label', 'Autofit (coarse fit, fine fit & plot)', 'CallBack', comauto);
    uimenu( submenu, 'Label', 'Plot component dipoles'   , 'CallBack', complot, 'separator', 'on');

author: PatternRecognition

repo: OpenBMI


text:
% this is a demo low-performance EPI sequence;
% it doesn't use ramp-samping and is only good for educational purposes.
%
seq=mr.Sequence();              % Create a new sequence object
fov=220e-3; Nx=64; Ny=64;       % Define FOV and resolution
thickness=3e-3;                 % slice thinckness
Nslices=3;

% Set system limits
lims = mr.opts('MaxGrad',32,'GradUnit','mT/m',...
               'MaxSlew',130,'SlewUnit','T/m/s', ...
               'rfRingdownTime', 30e-6, 'rfDeadTime', 100e-6);

% Create 90 degree slice selection pulse and gradient
[rf, gz] = mr.makeSincPulse(pi/2,'system',lims,'Duration',3e-3,...
    'SliceThickness',thickness,'apodization',0.5,'timeBwProduct',4);

% Define other gradients and ADC events
deltak=1/fov;
kWidth = Nx*deltak;
dwellTime = 4e-6; % I want it to be divisible by 2
readoutTime = Nx*dwellTime;
flatTime=ceil(readoutTime*1e5)*1e-5; % round-up to the gradient raster
gx = mr.makeTrapezoid('x',lims,'Amplitude',kWidth/readoutTime,'FlatTime',flatTime);
adc = mr.makeAdc(Nx,'Duration',readoutTime,'Delay',gx.riseTime+flatTime/2-(readoutTime-dwellTime)/2);

% Pre-phasing gradients
preTime=8e-4;
gxPre = mr.makeTrapezoid('x',lims,'Area',-gx.area/2,'Duration',preTime); % removed -deltak/2 to aligh the echo between the samples
gzReph = mr.makeTrapezoid('z',lims,'Area',-gz.area/2,'Duration',preTime);
gyPre = mr.makeTrapezoid('y',lims,'Area',-Ny/2*deltak,'Duration',preTime);

% Phase blip in shortest possible time
dur = ceil(2*sqrt(deltak/lims.maxSlew)/10e-6)*10e-6;
gy = mr.makeTrapezoid('y',lims,'Area',deltak,'Duration',dur);

% Define sequence blocks
% seq.addBlock(mr.makeDelay(1)); % older scanners like Trio may need this
                                 % dummy delay to keep up with timing
for s=1:Nslices
    rf.freqOffset=gz.amplitude*thickness*(s-1-(Nslices-1)/2);
    seq.addBlock(rf,gz);
    seq.addBlock(gxPre,gyPre,gzReph);
    for i=1:Ny
        seq.addBlock(gx,adc);           % Read one line of k-space
        seq.addBlock(gy);               % Phase blip
        gx.amplitude = -gx.amplitude;   % Reverse polarity of read gradient
    end
end

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

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

%% Plot sequence waveforms
seq.plot();

%% trajectory calculation
[ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation, t_refocusing] = seq.calculateKspacePP();

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

seq.write('epi.seq');   % Output sequence for scanner
% seq.sound(); % simulate the seq's tone

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function [bestEpsilon bestF1] = selectThreshold(yval, pval)
%SELECTTHRESHOLD Find the best threshold (epsilon) to use for selecting
%outliers
%   [bestEpsilon bestF1] = SELECTTHRESHOLD(yval, pval) finds the best
%   threshold to use for selecting outliers based on the results from a
%   validation set (pval) and the ground truth (yval).
%

bestEpsilon = 0;
bestF1 = 0;
F1 = 0;

stepsize = (max(pval) - min(pval)) / 1000;
for epsilon = min(pval):stepsize:max(pval)

    % ====================== YOUR CODE HERE ======================
    % Instructions: Compute the F1 score of choosing epsilon as the
    %               threshold and place the value in F1. The code at the
    %               end of the loop will compare the F1 score for this
    %               choice of epsilon and set it to be the best epsilon if
    %               it is better than the current choice of epsilon.
    %
    % Note: You can use predictions = (pval < epsilon) to get a binary vector
    %       of 0's and 1's of the outlier predictions

    cvPrediction = pval<epsilon;
    tp = sum((cvPrediction == 1) & (yval == 1));
    fp = sum((cvPrediction == 1) & (yval == 0));
    fn = sum((cvPrediction == 0) & (yval == 1));
    precision = tp/(tp+fp);
    recision = tp/(tp+fn);
    F1 = (2*precision*recision)/(precision+recision);

    % =============================================================

    if F1 > bestF1
       bestF1 = F1;
       bestEpsilon = epsilon;
    end
end

end

repo: MachineLearningEx


text:
function ff=ref_fac(f,W,c,d,p,q,permutation)
%REF_FAC  Reference factorization.
%
%  This function cannot handle multidimensional arrays.
%  Reshape to matrix BEFORE you call this.
%

% Output
ff=zeros(p*q*W,c*d);

% If d==1, it is not possible to further reduce the size of
% wk. Actually, some of the following code produces an
% error, because Matlab interprets an fft of a 1x q*W*p as a
% a row operation !
if d>1

  % Shape to match first fft pass.
  work=zeros(d,q*W*p);

  % This loop iterates over the number of truly different wk's.
  for ko=0:c-1

    % Permute and copy into work array.
    % Format is suited for fft.
    work(:)=f(permutation+ko,:);

    % Execute the fft and place transposed in ff.
    ff(:,1+ko*d:(ko+1)*d)=fft(work.',[],2);

  end;

else

  % Work arrays.
  work=zeros(p*q*W,1);

  for ko=0:c-1

    % Permute input.
    work(:)=f(permutation+ko,:);

    % Write the block.
    ff(:,ko+1)=work;
  end

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function p = mahdist(x, m, S, V)
% mahalanobis distance from x to m by S
% mahdist([x(:)';y(:)';z(:)'],mu,[],S)
%   See NORMPDF for argument description.
% [d, n] = size(x)
% [d, 1] = size(m)
% S = []
% V = cov
[d, n] = size(x);
if nargin == 1
  dx = x;
elseif isempty(m)
  dx = x;
else
  % m specified
  sz = size(m);
  if sz(1) ~= d
    error('rows(m) ~= rows(x)')
  end
  nm = sz(2);
  if nm == 1
    dx = x - repmat(m,1,n);
  elseif n == 1
    dx = repmat(x,1,nm) - m;
  elseif nm == n
    dx = x - m;
  else
    error('incompatible number of columns in x and m')
  end
end
if nargin < 3
  % unit variance
  p = col_sum(dx.*dx);
  return
end
have_inv = 0;
if nargin == 3
  % standard deviation given
  if d == 1
    dx = dx./S;
    p = dx.*dx;
    return;
  end
  if S(2,1) ~= 0
    error('S is not upper triangular')
  end
  if any(size(S) ~= [d d])
    error('S is not the right size')
  end
else
  if ischar(V)
    if strcmp(V,'inv')
      % inverse stddev given
      iS = S;
      have_inv = 1;
    else
      error('unknown directive')
    end
  elseif ischar(S)
    if strcmp(S,'inv')
      % inverse variance given
      if d == 1
	iS = sqrt(V);
      else
	iS = chol(V);
      end
      have_inv = 1;
    else
      error('unknown directive')
    end
  else
    % variance given
    if d == 1
      S = sqrt(V);
    else
      S = chol(V);
    end
  end
end
if have_inv
  if d == 1
    dx = iS .* dx;
  else
    dx = iS*dx;
  end
else
  if d == 1
    dx = dx./S;
  else
    dx = solve_tril(S',dx);
  end
end
p = sum(dx.*dx,1);

author: kristinbranson

repo: JAABA


text:
function [ds]=sap2000elementos14(filenam)
%%%%%%%%%%%%%%%%%%           SAP2000 v14           %%%%%%%%%%%%%%%%%%%%%
%programa registrado por MECAVIL ingenieria-ciencia
%bresler_lin@hotmail.com
%%%este programa enumera los elementos y su nudo inicial y final
%arc=fopen('quimica3.$2k','r');
%filenam
arc=fopen(filenam,'r');
arc;
encuentra =0;
while encuentra==0
    linea=fgets(arc);
    if not(isempty(findstr('CONNECTIVITY - FRAME',linea)));
        encuentra=1;
   end
end
 conecc=[0;0;0];
 linea=fgets(arc);
 f=1;
 encuentra=0;
 while encuentra==0;
     [s1 s2 s3 s4 ] = strread(linea,'%s %s %s %s');
      s1=char(s1);s2=char(s2);s3=char(s3);
      v=[s1,s2,s3];
 [e1 e2]=strread(v,'%s %f ', 'delimiter',' =')  ;
  conecc=[ conecc e2(:,1)];
     f=f+1;
     linea=fgets(arc);
     linea;
     if not (isempty(findstr('COORDINATE SYSTEMS',linea)))
         encuentra=1;
     end
 end
 fclose(arc);
   conecc(:,1)=[];
   conecc= conecc';
    conecc;
  [aw1 aw2]=size( conecc)

  ds=[0 0 0];
  for i=1:1:aw1;
      ds=[ds;conecc(i,:)];
  end
  ds(1,:)=[];
inifr=ds(:,2);
finfr=ds(:,3);
  save elemento.txt ds -ASCII  ;%%%%%LOPUSE por no corregir las inercias
    save conecctividad.txt ds -ASCII  ;
%save quimicainicio.txt inifr -ASCII  ;
%save quimicafin.txt finfr -ASCII  ;


text:
function [week,sec_of_week] = gps_time(julday)
% GPS_TIME   Conversion of Julian Day number to GPS week and
%	                Seconds of Week as reckoned from midnight between
%                  Saturday  and Sunday

% Written by Kai Borre
% January 7, 2016

    a = floor(julday+.5);
    b = a+1537;
    c = floor((b-122.1)/365.25);
    e = floor(365.25*c);
    f = floor((b-e)/30.6001);
    d = b-e-floor(30.6001*f)+rem(julday+.5,1);
    day_of_week = rem(floor(julday+.5),7);
    week = floor((julday-2444244.5)/7);
    % We add +1 as the GPS week starts at Saturday midnight
    sec_of_week = (rem(d,1)+day_of_week+1)*86400;
%%%%%%% end gps_time.m	%%%%%%%%%%%%%

author: yandld


text:
function SurfCov = surface_coverage2(Axis,Len,Vec,height,nl,ns)

% Computes surface coverage (number between 0 and 1) of points on cylinder
% surface defined by "Axis" and "Len". "Vec" are the vectors connecting
% points to the Axis and "height" are the heights of the points from
% the base of the cylinder

[U,W] = orthonormal_vectors(Axis);
Vec = Vec*[U W];
ang = atan2(Vec(:,2),Vec(:,1))+pi;
I = ceil(height/Len*nl);
I(I == 0) = 1;
I(I > nl) = nl;
J = ceil(ang/2/pi*ns);
J(J == 0) = 1;
K = [I J-1]*[1 nl]';
SurfCov = length(unique(K))/nl/ns;

author: InverseTampere

repo: TreeQSM


text:
classdef MTCMO < ALGORITHM
% <multi> <real/integer/label/binary/permutation> <constrained>
% Multitasking constrained multi-objective optimization

%------------------------------- Reference --------------------------------
% K. Qiao, K. Yu, B. Qu, J. Liang, H. Song, C. Yue, H. Lin, and K. C. Tan,
% Dynamic auxiliary task-based evolutionary multitasking for constrained
% multi-objective optimization, IEEE Transactions on Evolutionary
% Computation, 2022.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Kangjia Qiao

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population1 = Problem.Initialization();
            Fitness1    = CalFitness(Population1.objs,Population1.cons);

            Population2 = Problem.Initialization();
            Fitness2   = CalFitness(Population2.objs,Population2.cons);

            cons = [Population1.cons;Population2.cons];
            cons(cons<0) = 0;
            VAR0 = max(sum(cons,2));
            if VAR0 == 0
                VAR0 = 1;
            end
            X=0;

            %% Optimization
            while Algorithm.NotTerminated(Population1)
                %% Udate the epsilon value
                cp  = (-log(VAR0)-6)/log(1-0.5);
                VAR = VAR0*(1-X)^cp;
                %% Offspring generation
                MatingPool1 = TournamentSelection(2,Problem.N,Fitness1);
                Offspring1  = OperatorGAhalf(Problem,Population1(MatingPool1));

                MatingPool2 = TournamentSelection(2,Problem.N,Fitness2);
                Offspring2  = OperatorGAhalf(Problem,Population2(MatingPool2));
                %% Environmental selection
                [Population1,Fitness1] = Main_task_EnvironmentalSelection([Population1,Offspring1,Offspring2],Problem.N,true);
                [Population2,Fitness2] = Auxiliray_task_EnvironmentalSelection([Population2,Offspring2,Offspring1],Problem.N,VAR);
                X = X + 1/(Problem.maxFE/Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
moveRealtimeCircle()
%% The video associated with this demonstrator is at:
% https://youtu.be/ALBaCQzSeu0

% Copyright Mohammad SAFEEA, 2nd April 2018

repo: KST-Kuka-Sunrise-Toolbox


text:
function nScatt = diracnorm_scat(x,cascade)
	% Calculation of a normalized scattering representation
	% nScatt is an array of size 2N which stores a
	% piecewise constant signal giving
	% the 2^J normalized values || SJ [p]f||/||SJ[p] delta||
	% over intervals of lengths proportional to ||SJ[p] delta||^2 .
	% Normalized scattering transforms are defined in the paper
	% "Group Invariant Scattering", S. Mallat, http://arxiv.org/abs/1101.2286

	N=size(x,1); % N must be an integer multiplied by a very high power of 2
	%Compute the scattering transform of the signal x and output it in a table
	Sx=scat(x,cascade);

	%Prepare the dirac signal
	dirac = zeros(N,1);
	dirac(N/2) = 1;
	%Compute its scattering transform and output it in the form of a table

	Sdirac=scat(dirac, cascade);
	%%
	%fprintf('scat time\n');
	[Scatt order2]=reorder_scat(Sx);

	%%
	nSdirac=reorder_scat(Sdirac);

	%%
	%fprintf('reordering time\n');

	lsig=2*N/2^(length(Sx{2}.signal));

	for k=0:lsig:2*N-lsig

		normd=norm(nSdirac(k+1:k+lsig));
		normx=norm(Scatt(k+1:k+lsig))/normd;

		% the order is coded in the smallest digits of normx
		normx = 10^(-6) * floor(normx * 10^(6));
		Scatt(k+1:k+lsig) = normx + order2(k+1)* 10^(-7);

		nSdirac(k+1:k+lsig)=normd^2;

	end

	C=sum(nSdirac)/lsig;

	p=0;
	for k=0:lsig:2*N-lsig
		step=nSdirac(k+1) * 2 * N/C;

		nScatt(floor(p+1):round(p+step))=Scatt(k+1);
		p=p+step;
	end

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [Y,T]=derivs(x,dim)
%Syntax: [Y,T]=derivs(x,dim)
%____________________________
%
% The phase space reconstruction of a time series x whith the derivatives approach,
% in embedding dimension m.
%
% Y is the trajectory matrix in the reconstructed phase space.
% T is the phase space length.
% x is the time series.
% dim is the embedding dimension.
%
%
% Reference:
%
% Packard N H, Cruchfield J P, Farmer J D, Shaw R S (1980): Geometry from a Time
% Series. Physical Review Letters 45: 712-715
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 11 Mar 2001.

if nargin<1 | isempty(x)==1
   error('You should provide a time series.');
else
   % x must be a vector
   if min(size(x))>1
      error('Invalid time series.');
   end
   x=x(:);
   % N is the time series length
   N=length(x);
end

if nargin<2 | isempty(dim)==1
    dim=2;
else
    % dim must be scalar
    if sum(size(dim))>2
        error('dim must be scalar.');
    end
    % dim must be an integer
    if dim-round(dim)~=0
        error('dim must be an integer.');
    end
    % dim must be positive
    if dim<=0
        error('dim must be positive.');
    end
end

% Initialize the phase space
Y=zeros(N,dim);

% Phase space reconstruction with the derivatives approach
Y(:,1)=x;
for i=2:dim
    Y(:,i)=[diff(Y(:,i-1));0];
end

% Total points on phase space
T=N-(dim-1);

% Remove the meaningless points
Y(T+1:end,:)=[];

author: Sable

repo: mcbench-benchmarks


function res = test_nested_1

% This is the top level function.
% Define some variables.
a = 1;
b = 2;
x = 0;
y = 9;

% Display variables before call to fun1
disp('Before call to fun1:')
disp(['a = ' num2str(a)]);
disp(['b = ' num2str(b)]);
disp(['x = ' num2str(x)]);
disp(['y = ' num2str(y)]);

% Call fun1
x = fun1(x);

% Display variables after call to fun1
disp(' ');
disp('After call to fun1:')
disp(['a = ' num2str(a)]);
disp(['b = ' num2str(b)]);
disp(['x = ' num2str(x)]);
disp(['y = ' num2str(y)]);

   % Declare a nested function
   function res = fun1(y)
   res = y + a;
   a = a + 1;

   % Display variables inside call to fun1
   disp(' ');
   disp('Inside call to fun1:')
   disp(['a = ' num2str(a)]);
   disp(['b = ' num2str(b)]);
   disp(['x = ' num2str(x)]);
   disp(['y = ' num2str(y)]);

   end % function fun1

end % function test_nested_1