text:
function out = spm_run_factorial_design(job)
% SPM job execution function - factorial design specification
% Input:
% job    - harvested job data structure (see matlabbatch help)
% Output:
% out    - struct variable containing the path of the saved SPM.mat
%__________________________________________________________________________
% Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_run_factorial_design.m 7739 2019-12-02 14:00:18Z guillaume $

%--------------------------------------------------------------------------
% This function configures the design matrix (describing the general
% linear model), data specification, and other parameters necessary for
% the statistical analysis. These parameters are saved in a
% configuration file (SPM.mat) in the current directory, and are
% passed on to spm_spm.m (via the Estimate button) which estimates the
% design. Inference on these estimated parameters is then handled by the
% SPM results section.
%
% This function, that sets up the necessary SPM structures, has been
% largely cannibalised from spm_spm_ui.m which we have retained for
% developmental continuity.
%
% It has in common with spm_spm_ui.m its use of the I factor matrix,
% the H,C,B,G design matrix partitions, the sF,sCFI,CFIforms,sCC,CCforms,
% sGXcalc,sGloNorm,sGMsca option definition variables and use of the
% functions spm_DesMtx.m, spm_meanby.m and spm_get_vc.m.
%
% It departs from spm_spm_ui.m in that it does not use the design
% definition data structure D. Also, it uses the new SPM.factor field,
% which for the case of full factorial designs, is used to automatically
% generate contrasts testing for main effects and interactions.
%
% This function departs from spm_spm_ui.m in that it does not provide
% the same menu of design options (these were hardcoded in D). Instead it
% provides a number of options for simple designs (1) One-sample t-test,
% (2) Two-sample t-test, (3) Paired t-test and (4) Multiple regression.
% Two facilities are provided for specifying more complicated designs
% (5) Full-factorial and (6) Flexible-factorial. These should be able to
% specify all design options (and more) that were available in SPM2.
% For each of these design types one can additionally specify regressors
% using the `covariates' option.
%
% Options (5) and (6) differ in the efficiency (eg. number of key
% strokes/button presses) with which a given design can be specified. For
% example, one-way ANOVAs can be specified using either option, but (5) is
% usually more efficient.
%
% Full-factorial designs
% ______________________
%
% This option is best used when you wish to test for all main effects and
% interactions in one-way, two-way or three-way ANOVAs.
%
% Design specification proceeds in 2 stages. Firstly, by creating new
% factors and specifying the number of levels and name for each.
% Nonsphericity, ANOVA-by-factor (for PET data) and scaling options (for
% PET data) can also be specified at this stage. Secondly, scans are
% assigned separately to each cell. This accomodates unbalanced designs.
%
% For example, if you wish to test for a main effect in the population
% from which your subjects are drawn and have modelled that effect at the
% first level using K basis functions (eg. K=3 informed basis functions)
% you can use a one-way ANOVA with K-levels. Create a single factor with K
% levels and then assign the data to each cell eg. canonical, temporal
% derivative and dispersion derivative cells, where each cell is assigned
% scans from multiple subjects.
%
% SPM will automatically generate the contrasts necessary to test for all
% main effects and interactions.
%
% Flexible-factorial designs
% __________________________
%
% In this option the design matrix is created a block at a time. You can
% decide whether you wish each block to be a main effect or a (two-way)
% interaction.
%
% This option is best used for one-way, two-way or three-way ANOVAs but
% where you do not wish to test for all possible main effects and
% interactions. This is perhaps most useful for PET where there is usually
% not enough data to test for all possible effects. Or for 3-way ANOVAs
% where you do not wish to test for all of the two-way interactions. A
% typical example here would be a group-by-drug-by-task analysis where,
% perhaps, only (i) group-by-drug or (ii) group-by-task interactions are
% of interest. In this case it is only necessary to have two-blocks in the
% design matrix - one for each interaction. The three-way interaction can
% then be tested for using a contrast that computes the difference between
% (i) and (ii).
%
% Design specification then proceeds in 3 stages. Firstly, factors are
% created and names specified for each. Nonsphericity, ANOVA-by-factor and
% scaling options can also be specified at this stage.
%
% Secondly, a list of scans is produced along with a factor matrix, I.
% This is an nscan x 4 matrix of factor level indicators (see xX.I below).
% The first factor must be 'replication' but the other factors can be
% anything. Specification of I and the scan list can be achieved in one
% of two ways (a) the 'Specify All' option allows I to be typed in at the
% user interface or (more likely) loaded in from the matlab workspace.
% All of the scans are then selected in one go. (b) the 'Subjects' option
% allows you to enter scans a subject at a time. The corresponding
% experimental conditions (ie. levels of factors) are entered at the same
% time. SPM will then create the factor matrix I. This style of interface
% is similar to that available in SPM2.
%
% Thirdly, the design matrix is built up a block at a time. Each block
% can be a main effect or a (two-way) interaction.
%
%--------------------------------------------------------------------------
%
% Variables saved in the SPM stucture:
%
% xY.VY         - nScan x 1 struct array of memory mapped images
%                 (see spm_vol for definition of the map structure)
% xX            - structure describing design matrix
% xX.I          - nScan x 4 matrix of factor level indicators
%                 I(n,i) is the level of factor i corresponding to image n
% xX.sF         - 1x4 cellstr containing the names of the four factors
%                 xX.sF{i} is the name of factor i
% xX.X          - design matrix
% xX.xVi        - correlation constraints for non-spericity correction
% xX.iH         - vector of H partition (condition effects) indices,
%                 identifying columns of X corresponding to H
% xX.iC         - vector of C partition (covariates of interest) indices
% xX.iB         - vector of B partition (block effects) indices
% xX.iG         - vector of G partition (nuisance variables) indices
% xX.name       - p x 1 cellstr of effect names corresponding to columns
%                 of the design matrix
%
% xC            - structure array of covariate details
% xC(i).rc      - raw (as entered) i-th covariate
% xC(i).rcname  - name of this covariate (string)
% xC(i).c       - covariate as appears in design matrix (after any scaling,
%                 centering of interactions)
% xC(i).cname   - cellstr containing names for effects corresponding to
%                 columns of xC(i).c
% xC(i).iCC     - covariate centering option
% xC(i).iCFI    - covariate by factor interaction option
% xC(i).type    - covariate type: 1=interest, 2=nuisance, 3=global
% xC(i).cols    - columns of design matrix corresponding to xC(i).c
% xC(i).descrip - cellstr containing a description of the covariate
%
% xGX           - structure describing global options and values
% xGX.iGXcalc   - global calculation option used
% xGX.sGXcalc   - string describing global calculation used
% xGX.rg        - raw globals (before scaling and such like)
% xGX.iGMsca    - grand mean scaling option
% xGX.sGMsca    - string describing grand mean scaling
% xGX.GM        - value for grand mean (/proportional) scaling
% xGX.gSF       - global scaling factor (applied to xGX.rg)
% xGX.iGC       - global covariate centering option
% xGX.sGC       - string describing global covariate centering option
% xGX.gc        - center for global covariate
% xGX.iGloNorm  - Global normalisation option
% xGX.sGloNorm  - string describing global normalisation option
%
% xM            - structure describing masking options
% xM.T          - Threshold masking value (-Inf=>None, real=>absolute,
%                 complex=>proportional (i.e. times global))
% xM.TH         - nScan x 1 vector of analysis thresholds, one per image
% xM.I          - Implicit masking (0=>none, 1=>implicit zero/NaN mask)
% xM.VM         - struct array of explicit mask images
%                 (empty if no explicit masks)
% xM.xs         - structure describing masking options
%                 (format is same as for xsDes described below)
%
% xsDes         - structure of strings describing the design:
%                 Fieldnames are essentially topic strings (use "_"'s for
%                 spaces), and the field values should be strings or cellstr's
%                 of information regarding that topic. spm_DesRep.m
%                 uses this structure to produce a printed description
%                 of the design, displaying the fieldnames (with "_"'s
%                 converted to spaces) in bold as topics, with
%                 the corresponding text to the right
%
%--------------------------------------------------------------------------

%-Output directory
%--------------------------------------------------------------------------
cwd = pwd;
d   = spm_file(job.dir{1},'cpath');
if ~exist(d,'dir')
    sts = mkdir(d);
    if ~sts, error('Error creating output directory "%s".',d); end
end
cd(d);

%-Ask about overwriting files from previous analyses...
%--------------------------------------------------------------------------
if exist(fullfile(job.dir{1},'SPM.mat'),'file')
    str = { 'Current directory contains existing SPM file:',...
        'Continuing will overwrite existing file!'};
    if spm_input(str,1,'bd','stop|continue',[1,0],1,mfilename)
        fprintf('%-40s: %30s\n\n',...
            'Abort...   (existing SPM file)',spm('time'));
        return
    end
end

% If we've gotten to this point we're committed to overwriting files.
% Delete them so we don't get stuck in spm_spm
%--------------------------------------------------------------------------
files = {'^mask\..{3}$','^ResMS\..{3}$','^RPV\..{3}$',...
    '^beta_.{4}\..{3}$','^con_.{4}\..{3}$','^ResI_.{4}\..{3}$',...
    '^ess_.{4}\..{3}$', '^spm\w{1}_.{4}\..{3}$'};

for i=1:length(files)
    j = spm_select('List',pwd,files{i});
    for k=1:size(j,1)
        spm_unlink(deblank(j(k,:)));
    end
end

%-Option definitions
%==========================================================================

%-Generic factor names
%--------------------------------------------------------------------------
sF = {'sF1','sF2','sF3','sF4'};

%-Covariate by factor interaction options
sCFI = {'<none>';...                                        %-1
    'with sF1';'with sF2';'with sF3';'with sF4';...         %-2:5
    'with sF2 (within sF4)';'with sF3 (within sF4)'};       %-6,7

%-DesMtx argument components for covariate by factor interaction options
% (Used for CFI's Covariate Centering (CC), GMscale & Global normalisation)
%--------------------------------------------------------------------------
CFIforms = {'[]',   'C',    '{}';...                        %-1
    'I(:,1)',       'FxC',  '{sF{1}}';...                   %-2
    'I(:,2)',       'FxC',  '{sF{2}}';...                   %-3
    'I(:,3)',       'FxC',  '{sF{3}}';...                   %-4
    'I(:,4)',       'FxC',  '{sF{4}}';...                   %-5
    'I(:,[4,2])',   'FxC',  '{sF{4},sF{2}}';...             %-6
    'I(:,[4,3])',   'FxC',  '{sF{4},sF{3}}' };              %-7

%-Centre (mean correction) options for covariates & globals            (CC)
% (options 9-12 are for centering of global when using AnCova GloNorm) (GC)
%--------------------------------------------------------------------------
sCC = {'around overall mean';...                            %-1
    'around sF1 means';...                                  %-2
    'around sF2 means';...                                  %-3
    'around sF3 means';...                                  %-4
    'around sF4 means';...                                  %-5
    'around sF2 (within sF4) means';...                     %-6
    'around sF3 (within sF4) means';...                     %-7
    '<no centering>';...                                    %-8
    'around user specified value';...                       %-9
    '(as implied by AnCova)';...                            %-10
    'GM';...                                                %-11
    '(redundant: not doing AnCova)'}';                      %-12

%-DesMtx I forms for covariate centering options
%--------------------------------------------------------------------------
CCforms = {'ones(nScan,1)',CFIforms{2:end,1},''}';

%-Global calculation options                                       (GXcalc)
%--------------------------------------------------------------------------
sGXcalc  = {'omit';...                                      %-1
    'user specified';...                                    %-2
    'mean voxel value (within per image fullmean/8 mask)'}; %-3

%-Global normalization options                                    (GloNorm)
%--------------------------------------------------------------------------
sGloNorm = {'AnCova';...                                    %-1
    'AnCova by sF1';...                                     %-2
    'AnCova by sF2';...                                     %-3
    'AnCova by sF3';...                                     %-4
    'AnCova by sF4';...                                     %-5
    'AnCova by sF2 (within sF4)';...                        %-6
    'AnCova by sF3 (within sF4)';...                        %-7
    'proportional scaling';...                              %-8
    '<no global normalisation>'};                           %-9

%-Grand mean scaling options                                        (GMsca)
% (NB: Grand mean scaling by subject is redundent for proportional scaling)
%--------------------------------------------------------------------------
sGMsca = {'scaling of overall grand mean';...               %-1
    'scaling of sF1 grand means';...                        %-2
    'scaling of sF2 grand means';...                        %-3
    'scaling of sF3 grand means';...                        %-4
    'scaling of sF4 grand means';...                        %-5
    'scaling of sF2 (within sF4) grand means';...           %-6
    'scaling of sF3 (within sF4) grand means';...           %-7
    '(implicit in PropSca global normalisation)';...        %-8
    '<no grand Mean scaling>'   };                          %-9

%-Conditions of no interest defaults
%--------------------------------------------------------------------------
B      = [];
Bnames = {};
factor = [];

switch char(fieldnames(job.des))

    %-One sample t-test
    %======================================================================
    case 't1',

        DesName = 'One sample t-test';

        P = job.des.t1.scans;
        n = length(P);
        I = (1:n)';
        I = [I,ones(n,3)];

        [H,Hnames] = spm_DesMtx(I(:,2),'-','mean');

        factor(1).name     = 'Group';
        factor(1).levels   = 1;
        factor(1).variance = 0;
        factor(1).dept     = 0;

    %-Two-sample t-test
    %======================================================================
    case 't2',

        DesName = 'Two-sample t-test';

        P  = job.des.t2.scans1;
        n1 = length(job.des.t2.scans1);
        P  = [P; job.des.t2.scans2];
        n2 = length(job.des.t2.scans2);

        I  = [(1:n1),(1:n2)]';
        I  = [I,[ones(n1,1);2*ones(n2,1)]];
        I  = [I,ones(n1+n2,2)];

        [H,Hnames] = spm_DesMtx(I(:,2),'-','Group');

        % Names and levels
        factor(1).name     = 'Group';
        factor(1).levels   = 2;

        % Ancova options
        factor(1).gmsca    = job.des.t2.gmsca;
        factor(1).ancova   = job.des.t2.ancova;

        % Nonsphericity options
        factor(1).variance = job.des.t2.variance;
        factor(1).dept     = job.des.t2.dept;

        if any([n1 n2]==1) && factor(1).variance == 1
            warning('Imposing equal variance between groups for 1 vs N comparison.');
            factor(1).variance = 0;
        end

    %-Paired t-test
    %======================================================================
    case 'pt',

        DesName = 'Paired t-test';

        Npairs  = length(job.des.pt.pair);
        P       = [];
        for p   = 1:Npairs
            P   = [P;job.des.pt.pair(p).scans];
        end

        I       = ones(Npairs*2,1);
        I(:,2)  = kron([1:Npairs]',ones(2,1));
        I(:,3)  = kron(ones(Npairs,1),[1 2]');
        I(:,4)  = I(:,1);

        [B,Bnames] = spm_DesMtx(I(:,2),'-','Subject');
        [H,Hnames] = spm_DesMtx(I(:,3),'-','Condition');

        % Names and levels
        factor(1).name     = 'Subject';
        factor(1).levels   = Npairs;
        factor(2).name     = 'Condition';
        factor(2).levels   = 2;

        % Ancova options
        factor(1).gmsca    = 0;
        factor(1).ancova   = 0;
        factor(2).gmsca    = job.des.pt.gmsca;
        factor(2).ancova   = job.des.pt.ancova;

        % Nonsphericity options
        factor(1).variance = 0;
        factor(1).dept     = 0;
        factor(2).variance = 0;
        factor(2).dept     = 0;

    %-Multiple regression
    %======================================================================
    case 'mreg',

        DesName = 'Multiple regression';

        P = job.des.mreg.scans;
        n = length(P);
        I = (1:n)';
        I = [I,ones(n,3)];

        % Names and levels
        factor(1).name     = '';
        factor(1).levels   = 1;

        % Nonsphericity options
        factor(1).variance = 0;
        factor(1).dept     = 0;

        if job.des.mreg.incint==0
            H = []; Hnames = '';
        else
            [H,Hnames] = spm_DesMtx(I(:,2),'-','mean');
        end

        for i=1:length(job.des.mreg.mcov)
            job.cov(end+1).c   = job.des.mreg.mcov(i).c;
            job.cov(end).cname = job.des.mreg.mcov(i).cname;
            job.cov(end).iCC   = job.des.mreg.mcov(i).iCC;
            job.cov(end).iCFI  = 1;
        end

    %-ANOVA
    %======================================================================
    case 'anova',

        DesName = 'ANOVA';

        job.des.anova.fact.name = 'Groups';

        % Automatically number cells 1 to levels, so user doesn't have to
        levels = length(job.des.anova.icell);
        job.des.anova.fact.levels = levels;
        for i=1:levels
            job.des.anova.icell(i).levels = i;
        end
        [I,P,H,Hnames] = spm_design_factorial(job.des.anova);

        % Names and levels
        factor(1).name     = 'Groups';
        factor(1).levels   = levels;

        % Ancova options
        factor(1).gmsca    = job.des.anova.gmsca;
        factor(1).ancova   = job.des.anova.ancova;

        % Nonsphericity options
        factor(1).variance = job.des.anova.variance;
        factor(1).dept     = job.des.anova.dept;

    %-ANOVA: within-subject
    %======================================================================
    case 'anovaw',

        DesName = 'ANOVA - within subject';

        anovaw  = job.des.anovaw;
        anovaw.fac(1).name      = 'Subject';
        anovaw.fac(1).dept      = 0;
        anovaw.fac(1).variance  = 0;
        anovaw.fac(1).gmsca     = 0;
        anovaw.fac(1).ancova    = 0;

        anovaw.fac(2).name      = 'Groups';
        anovaw.fac(2).dept      = job.des.anovaw.dept;
        anovaw.fac(2).variance  = job.des.anovaw.variance;
        anovaw.fac(2).gmsca     = job.des.anovaw.gmsca;
        anovaw.fac(2).ancova    = job.des.anovaw.ancova;

        anovaw.fsuball.fsubject = anovaw.fsubject;

        % Main effect of subject and group
        anovaw.maininters{1}.fmain.fnum = 1;
        anovaw.maininters{2}.fmain.fnum = 2;

        [I,P,job.cov] = spm_design_within_subject(anovaw,job.cov);

        [H,Hnames,B,Bnames] = spm_design_flexible(anovaw,I);

        factor           = anovaw.fac;
        factor(1).levels = length(job.des.anovaw.fsubject);

    %-Full Factorial Design
    %======================================================================
    case 'fd',

        DesName = 'Full factorial';

        [I,P,H,Hnames] = spm_design_factorial(job.des.fd);

        Nfactors = length(job.des.fd.fact);
        for i=1:Nfactors
            % Names and levels
            factor(i).name     = job.des.fd.fact(i).name;
            factor(i).levels   = job.des.fd.fact(i).levels;

            % Ancova options
            factor(i).gmsca    = job.des.fd.fact(i).gmsca;
            factor(i).ancova   = job.des.fd.fact(i).ancova;

            % Nonsphericity options
            factor(i).variance = job.des.fd.fact(i).variance;
            factor(i).dept     = job.des.fd.fact(i).dept;
        end

    %-Flexible factorial design
    %======================================================================
    case 'fblock',

        DesName = 'Flexible factorial';

        if isfield(job.des.fblock.fsuball,'fsubject')
            % Data has been entered subject by subject
            nf = length(job.des.fblock.fac);
            [I,P,job.cov] = spm_design_within_subject(job.des.fblock,job.cov);
        else
            % Specify all scans and factor matrix
            I = job.des.fblock.fsuball.specall.imatrix;
            [ns,nI] = size(I);
            % Pad out factorial matrix to cover the four canonical factors
            if nI < 4
                warning('Padding factor matrix to have four columns.');
                I = [I, ones(ns,4-nI)];
            end
            % Get number of factors
            nf = length(job.des.fblock.fac);
            P  = job.des.fblock.fsuball.specall.scans;
        end

        if isempty(job.des.fblock.maininters)
            warning('No main effects or interactions have been specified.');
        end
        [H,Hnames,B,Bnames] = spm_design_flexible(job.des.fblock,I);

        for i=1:nf
            % Names and levels
            factor(i).name     = job.des.fblock.fac(i).name;
            factor(i).levels   = length(unique(I(:,i+1)));

            % Ancova options
            factor(i).gmsca    = job.des.fblock.fac(i).gmsca;
            factor(i).ancova   = job.des.fblock.fac(i).ancova;

            % Nonsphericity options
            factor(i).variance = job.des.fblock.fac(i).variance;
            factor(i).dept     = job.des.fblock.fac(i).dept;
        end

end

nScan = size(I,1);

%-Covariate partition(s): interest (C) & nuisance (G) excluding global
%==========================================================================
dstr = {'covariate','nuisance variable'};
C  = []; Cnames = {};                 %-Covariate DesMtx partitions & names
G  = []; Gnames = {};

xC = [];                              %-Struct array to hold raw covariates

%-Multiple covariates
%--------------------------------------------------------------------------
for m=1:numel(job.multi_cov)
    for n=1:numel(job.multi_cov(m).files)
        tmp   = load(job.multi_cov(m).files{n});
        names = {};
        if isstruct(tmp) % .mat
            if isfield(tmp,'R')
                R = tmp.R;
                if isfield(tmp,'names')
                    names = tmp.names;
                end
            else
                error(['Variable ''R'' not found in multiple ' ...
                    'covariates file ''%s''.'], job.multi_cov(m).files{n});
            end
        elseif isnumeric(tmp) % .txt
            R     = tmp;
            % read names from first line if commented?
        end
        for j=1:size(R,2)
            job.cov(end+1).c   = R(:,j);
            if isempty(names)
                job.cov(end).cname = sprintf('R%d%s',j);
            else
                job.cov(end).cname = names{j};
            end
            job.cov(end).iCFI  = job.multi_cov(m).iCFI;
            job.cov(end).iCC   = job.multi_cov(m).iCC;
        end
    end
end

%-Covariates
%--------------------------------------------------------------------------
nc = length(job.cov);                 %-Number of covariates
for i=1:nc

    c      = job.cov(i).c;
    cname  = job.cov(i).cname;
    rc     = c;                       %-Save covariate value
    rcname = cname;                   %-Save covariate name
    if job.cov(i).iCFI==1
        iCFI=1;
    else
        % SPMs internal factor numbers are 1 higher than specified in user
        % interface as, internally, the first factor is always `replication'
        iCFI = job.cov(i).iCFI+1;
    end
    switch job.cov(i).iCC
        case 1
            iCC = 1;
        case {2,3,4}
            iCC = job.cov(i).iCC + 1;
        otherwise
            iCC = job.cov(i).iCC + 3;
    end

    %-Centre within factor levels as appropriate
    if any(iCC == (1:7))
        c = c - spm_meanby(c,eval(CCforms{iCC}));
    end

    %-Do any interaction (only for single covariate vectors)
    %----------------------------------------------------------------------
    if iCFI > 1                       %-(NB:iCFI=1 if size(c,2)>1)
        tI        = [eval(CFIforms{iCFI,1}),c];
        tConst    = CFIforms{iCFI,2};
        tFnames   = [eval(CFIforms{iCFI,3}),{cname}];
        [c,cname] = spm_DesMtx(tI,tConst,tFnames);
    elseif size(c,2)>1                %-Design matrix block
        [null,cname] = spm_DesMtx(c,'X',cname);
    else
        cname     = {cname};
    end

    %-Store raw covariate details in xC struct for reference
    %-Pack c into appropriate DesMtx partition
    %----------------------------------------------------------------------
    %-Construct description string for covariate
    str = {sprintf('%s',rcname)};
    if size(rc,2)>1, str = {sprintf('%s (block of %d covariates)',...
            str{:},size(rc,2))}; end
    if iCC < 8, str=[str;{['used centered ',sCC{iCC}]}]; end
    if iCFI> 1, str=[str;{['fitted as interaction ',sCFI{iCFI}]}]; end

    typ = 1;
    tmp = struct(...
        'rc',   rc,    'rcname', rcname,...
        'c',    c,     'cname',  {cname},...
        'iCC',  iCC,   'iCFI',   iCFI,...
        'type', typ,...
        'cols', [1:size(c,2)] + size([H,C],2) + size([B,G],2)*min(typ-1,1),...
        'descrip', {str});
    if isempty(xC), xC = tmp; else xC = [xC,tmp]; end
    C     = [C,c];
    Cnames = [Cnames; cname];

end
clear c tI tConst tFnames

%==========================================================================
% - C O N F I G U R E   D E S I G N -
%==========================================================================

%-Images & image info: Map Y image files and check consistency of
% dimensions and orientation / voxel size
%==========================================================================
fprintf('%-40s: ','Mapping files')                                      %-#
VY    = spm_data_hdr_read(char(P));

%-Check compatibility of images
%--------------------------------------------------------------------------
spm_check_orientations(VY);

fprintf('%30s\n','...done')                                             %-#

%-Global values, scaling and global normalisation
%==========================================================================
%-Compute global values
%--------------------------------------------------------------------------
switch char(fieldnames(job.globalc))
    case 'g_omit',
        iGXcalc = 1;
    case 'g_user',
        iGXcalc = 2;
    case 'g_mean',
        iGXcalc = 3;
end

switch job.globalm.glonorm
    case 1,
        iGloNorm = 9;
    case 2,
        iGloNorm = 8;
    case 3,
        iGloNorm = 1;
end
if factor(1).levels > 1
    % Override if factor-specific ANCOVA has been specified
    for i=1:length(factor)
        if factor(i).ancova
            iGloNorm=i+2;
        end
    end
end

%-Analysis threshold mask
%--------------------------------------------------------------------------
%-Work out available options:
% -Inf=>None, real=>absolute, complex=>proportional, (i.e. times global)
M_T = -Inf;
switch char(fieldnames(job.masking.tm)),
    case 'tma',
        % Absolute
        M_T = job.masking.tm.tma.athresh;
    case 'tmr',
        % Relative
        M_T = job.masking.tm.tmr.rthresh*sqrt(-1);
        % Need to force calculation of globals
        if iGXcalc~=2, iGXcalc=3; end
    case 'tm_none'
        % None
        M_T = -Inf;
end

if iGXcalc==1 && (any(iGloNorm == [1:5 8]) || ...
        (factor(1).levels > 1 && any([factor.gmsca])))
    % Over-ride omission of global calculation if we need it
    disp(' ');
    disp('SPM needs estimates of global activity.');
    disp('But you have specified to omit this computation.');
    disp('SPM has overridden this omission and will automatically compute ');
    disp('globals as the mean value of within brain voxels.');
    disp(' ');
    iGXcalc = 3;
end
sGXcalc = sGXcalc{iGXcalc};

switch iGXcalc,
    case 1
        %-Don't compute => no GMsca (iGMsca==9) or GloNorm (iGloNorm==9)
        g = [];
    case 2
        %-User specified globals
        g = job.globalc.g_user.global_uval;
    case 3
        %-Compute as mean voxel value (within per image fullmean/8 mask)
        g = zeros(nScan,1);
        fprintf('%-40s: %30s','Calculating globals',' ')                %-#
        for i = 1:nScan
            str = sprintf('%3d/%-3d',i,nScan);
            fprintf('%s%30s',repmat(sprintf('\b'),1,30),str)            %-#
            g(i) = spm_global(VY(i)); % FIXME % for meshes
        end
        fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done')        %-#
    otherwise
        error('illegal iGXcalc')
end
rg = g;

fprintf('%-40s: ','Design configuration')                               %-#

%-Grand mean scaling options                                        (GMsca)
%--------------------------------------------------------------------------
if iGloNorm==8
    iGMsca=8;   %-grand mean scaling implicit in PropSca GloNorm
else
    switch char(fieldnames(job.globalm.gmsca))
        case 'gmsca_yes',
            iGMsca=1;
        case 'gmsca_no',
            iGMsca=9;
    end
    if factor(1).levels > 1
        % Over-ride if factor-specific scaling has been specified
        for i=1:numel(factor)
            if factor(i).gmsca
                iGMsca=i+2;
            end
        end
    end
end

%-Value for PropSca / GMsca                                            (GM)
%--------------------------------------------------------------------------
switch iGMsca,
    case 9                                %-Not scaling (GMsca or PropSca)
        GM = 0;                           %-Set GM to zero when not scaling
    case 1                                %-Ask user value of GM
        GM = job.globalm.gmsca.gmsca_yes.gmscv;
    otherwise
        if iGloNorm==8
            switch char(fieldnames(job.globalm.gmsca))
                case 'gmsca_yes',
                    % Proportionally scale to this value
                    GM = job.globalm.gmsca.gmsca_yes.gmscv;
                case 'gmsca_no',
                    GM = 50;
            end
        else
            % Grand mean scaling by factor eg. scans are scaled so that the
            % mean global value over each level of the factor is set to GM
            GM=50;
        end
end

%-If GM is zero then don't GMsca! or PropSca GloNorm
if GM==0,
    iGMsca=9;
    if iGloNorm==8,
        iGloNorm=9;
    end
end

%-Sort out description strings for GloNorm and GMsca
%--------------------------------------------------------------------------
sGloNorm = sGloNorm{iGloNorm};
sGMsca   = sGMsca{iGMsca};
if iGloNorm==8
    sGloNorm = sprintf('%s to %-4g',sGloNorm,GM);
elseif iGMsca<8
    sGMsca   = sprintf('%s to %-4g',sGMsca,GM);
end

%-Scaling: compute global scaling factors gSF required to implement
% proportional scaling global normalisation (PropSca) or grand mean
% scaling (GMsca), as specified by iGMsca (& iGloNorm)
%--------------------------------------------------------------------------
switch iGMsca,
    case 8
        %-Proportional scaling global normalisation
        if iGloNorm~=8, error('iGloNorm-iGMsca(8) mismatch for PropSca'), end
        gSF    = GM./g;
        g      = GM*ones(nScan,1);
    case {1,2,3,4,5,6,7}
        %-Grand mean scaling according to iGMsca
        if iGXcalc==1, error('Global calculation option is not appropriate.'), end
        gSF    = GM./spm_meanby(g,eval(CCforms{iGMsca}));
        g      = g.*gSF;
    case 9
        %-No grand mean scaling
        gSF    = ones(nScan,1);
    otherwise
        error('illegal iGMsca')
end

%-Apply gSF to memory-mapped scalefactors to implement scaling
%--------------------------------------------------------------------------
for i = 1:nScan
    VY(i).pinfo(1:2,:) = VY(i).pinfo(1:2,:)*gSF(i); % FIXME % for meshes
end

%-Global centering (for AnCova GloNorm)                                (GC)
%-If not doing AnCova then GC is irrelevant
%--------------------------------------------------------------------------
if ~any(iGloNorm == [1:7])
    iGC = 12;
    gc  = [];
else
    iGC = 10;
    gc = 0;
end

%-AnCova: Construct global nuisance covariates partition (if AnCova)
%--------------------------------------------------------------------------
if any(iGloNorm == [1:7])

    %-Centre global covariate as requested
    %----------------------------------------------------------------------
    switch iGC, case {1,2,3,4,5,6,7}    %-Standard sCC options
        gc = spm_meanby(g,eval(CCforms{iGC}));
        case 8                  %-No centering
            gc = 0;
        case 9                  %-User specified centre
            %-gc set above
        case 10                 %-As implied by AnCova option
            gc = spm_meanby(g,eval(CCforms{iGloNorm}));
        case 11                 %-Around GM
            gc = GM;
        otherwise               %-unknown iGC
            error('unexpected iGC value')
    end

    %-AnCova - add scaled centred global to DesMtx `G' partition
    %----------------------------------------------------------------------
    rcname     = 'global';
    tI         = [eval(CFIforms{iGloNorm,1}),g - gc];
    tConst     = CFIforms{iGloNorm,2};
    tFnames    = [eval(CFIforms{iGloNorm,3}),{rcname}];
    [f,gnames]  = spm_DesMtx(tI,tConst,tFnames);
    clear tI tConst tFnames

    %-Save GX info in xC struct for reference
    %----------------------------------------------------------------------
    str     = {sprintf('%s: %s',dstr{2},rcname)};
    if any(iGMsca==[1:7]), str=[str;{['(after ',sGMsca,')']}]; end
    if iGC ~= 8, str=[str;{['used centered ',sCC{iGC}]}]; end
    if iGloNorm > 1
        str=[str;{['fitted as interaction ',sCFI{iGloNorm}]}];
    end
    tmp  = struct(  'rc',rg.*gSF,       'rcname',rcname,...
        'c',f,          'cname' ,{gnames},...
        'iCC',iGC,      'iCFI'  ,iGloNorm,...
        'type',         3,...
        'cols',[1:size(f,2)] + size([H C B G],2),...
        'descrip',      {str}       );

    G = [G,f]; Gnames = [Gnames; gnames];
    if isempty(xC), xC = tmp; else xC = [xC,tmp]; end

elseif iGloNorm==8 || iGXcalc>1

    %-Globals calculated, but not AnCova: Make a note of globals
    %----------------------------------------------------------------------
    if iGloNorm==8
        str = { 'global values: (used for proportional scaling)';...
            '("raw" unscaled globals shown)'};
    elseif isfinite(M_T) && ~isreal(M_T)
        str = { 'global values: (used to compute analysis threshold)'};
    else
        str = { 'global values: (computed but not used)'};
    end

    rcname ='global';
    tmp     = struct('rc',rg,    'rcname',rcname,...
        'c',{[]},   'cname' ,{{}},...
        'iCC',0,    'iCFI'  ,0,...
        'type',     3,...
        'cols',     {[]},...
        'descrip',  {str}           );

    if isempty(xC), xC = tmp; else xC = [xC,tmp]; end
end

%-Save info on global calculation in xGX structure
%--------------------------------------------------------------------------
xGX = struct(...
    'iGXcalc', iGXcalc,  'sGXcalc', sGXcalc,  'rg',rg,...
    'iGMsca',  iGMsca,   'sGMsca',  sGMsca,   'GM',GM,    'gSF',gSF,...
    'iGC',     iGC,      'sGC',     sCC{iGC}, 'gc',gc,...
    'iGloNorm',iGloNorm, 'sGloNorm',sGloNorm);

%-Make a description string
%--------------------------------------------------------------------------
if isinf(M_T)
    xsM.Analysis_threshold = 'None (-Inf)';
elseif isreal(M_T)
    xsM.Analysis_threshold = sprintf('images thresholded at %6g',M_T);
else
    xsM.Analysis_threshold = sprintf(['images thresholded at %6g ',...
        'times global'],imag(M_T));
end

%-Construct masking information structure and compute actual analysis
% threshold using scaled globals (rg.*gSF)
%--------------------------------------------------------------------------
if isreal(M_T),
    M_TH = M_T  * ones(nScan,1);    %-NB: -Inf is real
else
    M_TH = imag(M_T) * (rg.*gSF);
end

%-Implicit masking: Ignore zero voxels in low data-types?
%--------------------------------------------------------------------------
% (Implicit mask is NaN in higher data-types.)
if ~spm_type(VY(1).dt(1),'nanrep')
    M_I = job.masking.im;  % Implicit mask ?
    if M_I
        xsM.Implicit_masking = 'Yes: zero''s treated as missing';
    else
        xsM.Implicit_masking = 'No';
    end
else
    M_I = 1;
    xsM.Implicit_masking = 'Yes: NaN''s treated as missing';
end

%-Explicit masking
%--------------------------------------------------------------------------
if isempty(job.masking.em{:})
    VM = [];
    xsM.Explicit_masking = 'No';
else
    VM = spm_data_hdr_read(char(job.masking.em));
    xsM.Explicit_masking = 'Yes';
end

xM     = struct('T',M_T, 'TH',M_TH, 'I',M_I, 'VM',{VM}, 'xs',xsM);

%-Construct full design matrix (X), parameter names and structure (xX)
%==========================================================================
X      = [H C B G];
tmp    = cumsum([size(H,2), size(C,2), size(B,2), size(G,2)]);
xX     = struct(...
    'X',        X,...
    'iH',       [1:size(H,2)],...
    'iC',       [1:size(C,2)] + tmp(1),...
    'iB',       [1:size(B,2)] + tmp(2),...
    'iG',       [1:size(G,2)] + tmp(3),...
    'name',     {[Hnames; Cnames; Bnames; Gnames]},...
    'I',        I,...
    'sF',       {sF});

%-Design description (an nx2 cellstr) - for saving and display
%==========================================================================
tmp = {sprintf('%d condition, +%d covariate, +%d block, +%d nuisance',...
    size(H,2),size(C,2),size(B,2),size(G,2));...
    sprintf('%d total, having %d degrees of freedom',...
    size(X,2),rank(X));...
    sprintf('leaving %d degrees of freedom from %d images',...
    size(X,1)-rank(X),size(X,1))};
xsDes = struct('Design',    {DesName},...
    'Global_calculation',   {sGXcalc},...
    'Grand_mean_scaling',   {sGMsca},...
    'Global_normalisation', {sGloNorm},...
    'Parameters',           {tmp});

fprintf('%30s\n','...done')                                             %-#

%-Generate error covariance components (non-sphericity)
%==========================================================================
Vi          = spm_get_vc(I, factor);

%-Assemble SPM structure
%==========================================================================
SPM.xY.P    = P;            % filenames
SPM.xY.VY   = VY;           % mapped data
SPM.nscan   = size(xX.X,1); % scan number
SPM.xX      = xX;           % design structure
SPM.xC      = xC;           % covariate structure
SPM.xGX     = xGX;          % global structure
SPM.xM      = xM;           % mask structure
SPM.xsDes   = xsDes;        % description
SPM.xVi.I   = I;            % factor matrix
if numel(Vi) == 1
    SPM.xVi.V  = Vi{1};     % non-sphericity matrix
else
    SPM.xVi.Vi = Vi;        % non-sphericity variance components
end

%-Automatic contrast generation for 'Full factorial'
%--------------------------------------------------------------------------
if strcmp(char(fieldnames(job.des)),'fd') && job.des.fd.contrasts
    SPM.factor = factor;
end

%-Save SPM.mat and set output argument
%--------------------------------------------------------------------------
fprintf('%-40s: ','Saving SPM configuration')                           %-#
fmt = spm_get_defaults('mat.format');
s = whos('SPM');
if s.bytes > 2147483647, fmt = '-v7.3'; end
save('SPM.mat', 'SPM', fmt);
fprintf('%30s\n','...SPM.mat saved')                                    %-#

out.spmmat{1} = fullfile(pwd, 'SPM.mat');

%-Display Design report
%==========================================================================
if ~spm('CmdLine') && ~isempty(spm_figure('FindWin','Graphics'))
    fprintf('%-40s: ','Design reporting')                               %-#
    fname     = cat(1,{SPM.xY.VY.fname}');
    spm_DesRep('DesMtx',SPM.xX,fname,SPM.xsDes)
    fprintf('%30s\n','...done')                                         %-#
end

fprintf('%-40s: %30s\n','Completed',spm('time'))                        %-#

%-Change back directory
%--------------------------------------------------------------------------
cd(cwd);

author: spm


text:
function x = cifft2(xf)

if isa(xf, 'gpuArray')
    x = real(ifft2(ifftshift(ifftshift(xf, 1), 2)));
else
    x = ifft2(ifftshift(ifftshift(xf, 1), 2), 'symmetric');
end

author: martin-danelljan

repo: ECO


text:
function u = bc_test01 ( x_num, x, t, u )

%*****************************************************************************80
%
%% BC_TEST01 sets the boundary conditions for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X_NUM, the number of nodes.
%
%    Input, real  X(X_NUM,1), the coordinates of the points.
%
%    Input, real T, the current time.
%
%    Input, real U(X_NUM,1), the solution at time T.
%
%    Output, real U(X_NUM,1), the solution at time T, with
%    boundary conditions enforced.
%
  u(1,1)     = 90.0;
  u(x_num,1) = 70.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% p(y_n) = N(y_n | C*x_n, R)
% p(x_n | x_(n-1)) = N(x_n | A*x_(n-1), Q)
%
% [x, Covx, logp] = kalman_filter_step(x, Covx, y, A, Q, C, R, 1)
%
% [x, Covx, logp] = kalman_filter_step(x, Covx, y, A, Q, C, inv(chol(R,'lower')), 2)

% Last modified 2011-10-19
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@aalto.fi)

function [x, Covx, logp] = kalman_filter_step(x, Covx, y, A, Q, C, R, method)

% Prediction step
x = A*x;
Covx = A*Covx*A' + Q;

% Update step
switch method
 case 1

  % Standard update step (if length(y) <= length(x))
  y = y - C*x;
  S = C*Covx*C' + R;
  L = chol(S, 'lower');
  K = linsolve_tril(L, C*Covx); %(Covx * C') / S;
  v = linsolve_tril(L, y);
  if nargout >= 3
    logp = gaussian_logpdf(v'*v, ...
                           0, ...
                           0, ...
                           logdet_chol(L), ...
                           length(y));
  end
  x = x + K'*v;
  Covx = Covx - K'*K; % K*C*Covx;

 case 2
  % Update step in lower-dimensional space of X

  invL_R = R;
  Lx = chol(Covx, 'lower');
  Z = invL_R * C * Lx;
  L = chol(speye(size(Lx)) + Z'*Z, 'lower');
  v_y = y - C*x;
  v_x = C'*(invL_R'*(invL_R*y)) + linsolve_lchol(Lx,x);
  Z = linsolve_tril(L, Lx');
  Covx = Z'*Z;
  x = Covx * v_x;

  if nargout >= 3
    v = invL_R' * (invL_R*v_y); % R \ (y - C*x)
    logp = gaussian_logpdf(v'*v_y - v'*C*Covx*C'*v, ...
                           0, ...
                           0, ...
                           logdet_chol(L) - logdet_chol(invL_R), ...
                           length(y));
  end

% $$$   invR = R;
% $$$   S = Covx * C';
% $$$   K = Covx + S*invR*S';
% $$$   L = chol(K, 'lower');
% $$$   Z = L \ Covx;
% $$$   if nargout >= 3
% $$$     % Compute marginal likelihood term
% $$$   end
% $$$   x = Z' * linsolve_tril(L, S*(invR*y) + x);
% $$$   Covx = Z'*Z;

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [cscale, rscale] = gmscale(A, iprint, scltol)
% Geometric-Mean Scaling finds the scale values for the
% `m x n` sparse matrix `A`.
%
% USAGE:
%
%    [cscale, rscale] = gmscale(A, iprint, scltol)
%
% INPUTS:
%    A(i, j):           contains entries of `A`.
%    iprint:            > 0 requests messages to the screen (0 means no output).
%    scltol:            should be in the range (0.0, 1.0).
%                       Typically `scltol` = 0.9.  A bigger value like 0.99 asks
%                       `gmscale` to work a little harder (more passes).
%
% OUTPUTS:
%    cscale, rscale:    column vectors of column and row scales such that
%                       `R` (inverse) `A` `C` (inverse) should have entries near 1.0,
%                       where `R= diag(rscale)`, `C = diag(cscale)`.
%
% An iterative procedure based on geometric means is used,
% following a routine written by Robert Fourer, 1979.
% Several passes are made through the columns and rows of `A`.
% The main steps are:
%
%   1. Compute :math:`aratio = max_j (max_i Aij / min_i Aij)`.
%   2. Divide each row `i` by :math:`sqrt( max_j Aij * min_j Aij)`.
%   3. Divide each column `j` by :math:`sqrt( max_i Aij * min_i Aij)`.
%   4. Compute `sratio` as in Step 1.
%   5. If :math:`sratio < aratio * scltol`,
%      set :math:`aratio = sratio` and repeat from Step 2.
%
% To dampen the effect of very small elements, on each pass,
% a new row or column scale will not be smaller than sqrt(damp)
% times the largest (scaled) element in that row or column.
%
% Use of the scales:
% To apply the scales to a linear program,
% :math:`min c^T x` st :math:`A x = b`, :math:`l \leq x \leq u`,
% we need to define "barred" quantities by the following relations:
% `A = R Abar C`, `b = R bbar`, `C cbar = c`,
% `C l = lbar`, `C u = ubar`, `C x = xbar`.
%
% This gives the scaled problem
% :math:`min\ cbar^T xbar` st :math:`Abar\ xbar = bbar`, :math:`lbar \leq xbar \leq ubar`.
%
% .. Author: - Michael Saunders, Systems Optimization Laboratory, Stanford University.
% ..
%    07 Jun 1996: First f77 version, based on MINOS 5.5 routine m2scal.
%    24 Apr 1998: Added final pass to make column norms = 1.
%    18 Nov 1999: Fixed up documentation.
%    26 Mar 2006: (Leo Tenenblat) First Matlab version based on Fortran version.
%    21 Mar 2008: (MAS) Inner loops j = 1:n optimized.
%    09 Apr 2008: (MAS) All loops replaced by sparse-matrix operations.
%                 We can't find the biggest and smallest Aij
%                 on each scaling pass, so no longer print them.
%    24 Apr 2008: (MAS, Kaustuv) Allow for empty rows and columns.
%    13 Nov 2009: gmscal.m renamed gmscale.m.

  if iprint > 0
    fprintf('\ngmscale: Geometric-Mean scaling of matrix')
    fprintf('\n-------\n                 Max col ratio')
  end

  [m,n]   = size(A);
  A       = abs(A);    % Work with |Aij|
  maxpass = 10;
  aratio  = 1e+50;
  damp    = 1e-4;
  small   = 1e-8;
  rscale  = ones(m,1);
  cscale  = ones(n,1);

%---------------------------------------------------------------
% Main loop.
%---------------------------------------------------------------
  for npass = 0:maxpass

    % Find the largest column ratio.
    % Also set new column scales (except on pass 0).

    rscale(rscale==0) = 1;
    Rinv    = diag(sparse(1./rscale));
    SA      = Rinv*A;
    [I,J,V] = find(SA);
    invSA   = sparse(I,J,1./V,m,n);
    cmax    = full(max(   SA))';   % column vector
    cmin    = full(max(invSA))';   % column vector
    cmin    = 1./(cmin + eps);
    sratio  = max( cmax./cmin );   % Max col ratio
    if npass > 0
      cscale = sqrt( max(cmin, damp*cmax) .* cmax );
    end

    if iprint > 0
      fprintf('\n  After %2g %19.2f', npass, sratio)
    end

    if npass >= 2 && sratio >= aratio*scltol, break; end
    if npass == maxpass, break; end
    aratio  = sratio;

    % Set new row scales for the next pass.

    cscale(cscale==0) = 1;
    Cinv    = diag(sparse(1./cscale));
    SA      = A*Cinv;                  % Scaled A
    [I,J,V] = find(SA);
    invSA   = sparse(I,J,1./V,m,n);
    rmax    = full(max(   SA,[],2));   % column vector
    rmin    = full(max(invSA,[],2));   % column vector
    rmin    = 1./(rmin + eps);
    rscale  = sqrt( max(rmin, damp*rmax) .* rmax );
  end
%---------------------------------------------------------------
% End of main loop.
%---------------------------------------------------------------

% Reset column scales so the biggest element
% in each scaled column will be 1.
% Again, allow for empty rows and columns.

  rscale(rscale==0) = 1;
  Rinv    = diag(sparse(1./rscale));
  SA      = Rinv*A;
  [I,J,V] = find(SA);
  cscale  = full(max(SA))';   % column vector
  cscale(cscale==0) = 1;

% Find the min and max scales.

  if iprint>0
    [rmin,imin] = min(rscale);
    [rmax,imax] = max(rscale);
    [cmin,jmin] = min(cscale);
    [cmax,jmax] = max(cscale);

    fprintf('\n\n  Min scale               Max scale')
    fprintf('\n  Row %6g %9.1e    Row %6g %9.1e'  , imin, rmin, imax, rmax)
    fprintf('\n  Col %6g %9.1e    Col %6g %9.1e\n', jmin, cmin, jmax, cmax)
  end

% end of gmscale

author: opencobra

repo: cobratoolbox


text:
function [imgStats] = computeImageFeatures(inputImages, varargin)
	% [imgStats] = computeImageFeatures(inputImages, varargin)
	%
	% Filters large and small objects in an set of images, returns filtered matricies along with vector with decisions and sizes.
	%
	% Biafra Ahanonu
	% 2013.10.31
	% based on SpikeE code
	%
	% inputs
	%   inputImages - [x y nSignals]
	%
	% outputs
	%   imgStats -
	%
	% options
	%   minNumPixels
	%   maxNumPixels
	%   thresholdImages

	% changelog
		% updated: 2013.11.08 [09:24:12] removeSmallICs now calls a filterImages, name-change due to alteration in function, can slowly replace in codes
		% 2017.01.14 [20:06:04] - support switched from [nSignals x y] to [x y nSignals]
		% 2019.07.17 [00:29:16] - Added support for sparse input images (mainly ndSparse format).
		% 2019.09.10 [20:51:00] - Converted to parfor and removed unpacking options, bad practice.
		% 2019.10.02 [21:25:36] - Updated addedFeatures support for parfor and allowed this feature to be properly accessed by users.
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2022.04.13 [02:28:15] - Make default thresholding fast thresholding.
		% 2022.07.20 [14:42:00] - Improved annotation of code and added options.fastThresholding. Also by default plots are not made. Refactored extra features code and finding centroid to avoid unnecessary function calls, speeding up code.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% get options
	options.minNumPixels = 25;
	options.maxNumPixels = 600;
	options.makePlots = 0;
	options.waitbarOn = 1;
	options.thresholdImages = 1;
	options.threshold = 0.5;
	options.valid = [];
	% options.featureList = {'Eccentricity','EquivDiameter','Area','Orientation','Perimeter','Solidity',};
	options.featureList = {'Eccentricity','EquivDiameter','Area','Orientation','Perimeter'};
	% Whether to calculate non-regionprops
	options.addedFeatures = 0;
	% Input images for add features
	options.addedFeaturesInputImages = [];
	options.runRegionprops = 1;
	%
	options.xCoords = [];
	options.yCoords = [];
	%
	options.parforAltSwitch = 0;
	% Binary: 1 = fast thresholding (vectorized), 0 = normal thresholding
	options.fastThresholding = 1;
	% image filter: none, median,
	options.imageFilter = 'none';

	% Binary: 1 = whether to display info on command line.
	options.displayInfo = 1;

	options = getOptions(options,varargin);
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	parforAltSwitch = options.parforAltSwitch;
	% switch nargin
	% 	case 2
	% 		parforAltSwitch = 1;
	% 	case 1
	% 		% do nothing
	% 	otherwise
	% 		% do nothing
	% end

	nImages = size(inputImages,3);

	% reverseStr = '';
	% decide whether to threshold images
	if options.thresholdImages==1
		if options.addedFeatures==1&&isempty(options.addedFeaturesInputImages)
			options.addedFeaturesInputImages = inputImages;
		end
		inputImages = thresholdImages(inputImages,'waitbarOn',1,'binary',1,'threshold',options.threshold,'fastThresholding',options.fastThresholding,'removeUnconnected',1,'imageFilter',options.imageFilter);
	end

	% Only implement in Matlab 2017a and above
	if ~verLessThan('matlab', '9.2')
		D = parallel.pool.DataQueue;
		afterEach(D, @nUpdateParforProgress);
		p = 1;
		% N = nSignals;
	end
	nInterval = round(nImages/20);

	options_runRegionprops = options.runRegionprops;
	options_featureList = options.featureList;
	options_waitbarOn = options.waitbarOn;

	% loop over images and get their stats
	if options.displayInfo==1
		disp('Computing image features...')
	end
	regionStat = cell([nImages 1]);
	if options_runRegionprops==1
		% ticBytes(gcp)
		parfor imageNo = 1:nImages
			% iImage = squeeze(inputImages(:,:,imageNo));
			% imagesc(iImage)
			% imgStats.imageSizes(imageNo) = sum(iImage(:)>0);
			thisFilt = inputImages(:,:,imageNo);
			if issparse(thisFilt)
				thisFilt = full(thisFilt);
			end
			regionStat{imageNo} = regionprops(thisFilt, options_featureList);
			if ~verLessThan('matlab', '9.2')
				send(D, imageNo); % Update
			end
		end
		% tocBytes(gcp)
		% whos

		% Add region states to the general pool
		if ~verLessThan('matlab', '9.2'); p=1;end
		if options.displayInfo==1
			disp('Adding regionprops features to output...')
		end
		for imageNo = 1:nImages
			for ifeature = options_featureList
				% regionStat = regionprops(iImage, ifeature{1});
				try
					% eval(['imgStats.' ifeature{1} '(imageNo) = regionStat.' ifeature{1} ';']);
					imgStats.(ifeature{1})(imageNo) = regionStat{imageNo}.(ifeature{1});
				catch
					% eval(['imgStats.' ifeature{1} '(imageNo) = NaN;']);
					imgStats.(ifeature{1})(imageNo) = NaN;
				end
			end
			if ~verLessThan('matlab', '9.2')
				send(D, imageNo); % Update
			end
		end
	end

	if ~verLessThan('matlab', '9.2'); p=1;end

	if options.addedFeatures==1
		% get the centroids and other info for movie
		if isempty(options.xCoords)
			[xCoords, yCoords] = findCentroid(inputImages,'waitbarOn',options.waitbarOn,'runImageThreshold',0);
		else
			xCoords = options.xCoords;
			yCoords = options.yCoords;
		end

		if parforAltSwitch==1
			imgKurtosis = NaN([1 nImages]);
			imgSkewness = NaN([1 nImages]);
			if options.displayInfo==1
				disp('Computing alternative image features with parfor...')
			end
			addedFeaturesInputImages = options.addedFeaturesInputImages;
			parfor imageNo = 1:nImages
				thisFilt = addedFeaturesInputImages(:,:,imageNo);
				if issparse(thisFilt)
					thisFilt = full(thisFilt);
				end
				if isnan(xCoords(imageNo))==0
					t1 = getObjCutMovie(thisFilt,thisFilt,'cropSize',30,'createMontage',0,'crossHairsOn',0,'addPadding',1,'waitbarOn',0,'xCoords',xCoords(imageNo),'yCoords',yCoords(imageNo));
					t1 = t1{1};
					imgKurtosis(imageNo) = double(kurtosis(t1(:)));
					imgSkewness(imageNo) = double(skewness(t1(:)));
				else
					imgKurtosis(imageNo) = NaN;
					imgSkewness(imageNo) = NaN;
				end
				if ~verLessThan('matlab', '9.2')
					send(D, imageNo); % Update
				end
			end

			imgStats.imgKurtosis = imgKurtosis;
			imgStats.imgSkewness = imgSkewness;
		else
			if options.displayInfo==1
				disp('Computing alternative image features...')
			end
			for imageNo = 1:nImages
				% iImage2 = squeeze(options.addedFeaturesInputImages(:,:,imageNo));
				% figure(11);imagesc(iImage2);title(num2str(imageNo))
				% [imageNo xCoords(imageNo) yCoords(imageNo)]
				thisFilt = options.addedFeaturesInputImages(:,:,imageNo);
				if issparse(thisFilt)
					thisFilt = full(thisFilt);
				end
				if isnan(xCoords(imageNo))==0
					% t1=getObjCutMovie(options.addedFeaturesInputImages(:,:,imageNo),options.addedFeaturesInputImages(:,:,imageNo),'cropSize',30,'createMontage',0,'crossHairsOn',0,'addPadding',1,'waitbarOn',0,'xCoords',xCoords(imageNo),'yCoords',yCoords(imageNo));
					t1 = getObjCutMovie(thisFilt,thisFilt,'cropSize',30,'createMontage',0,'crossHairsOn',0,'addPadding',1,'waitbarOn',0,'xCoords',xCoords(imageNo),'yCoords',yCoords(imageNo));
					t1 = t1{1};
					imgStats.imgKurtosis(imageNo) = double(kurtosis(t1(:)));
					imgStats.imgSkewness(imageNo) = double(skewness(t1(:)));
				else
					imgStats.imgKurtosis(imageNo) = NaN;
					imgStats.imgSkewness(imageNo) = NaN;
				end
				% imgStats.imgKurtosis(imageNo) = kurtosis(iImage(:));
				% imgStats.imgSkewness(imageNo) = skewness(iImage(:));
				if ~verLessThan('matlab', '9.2')
					send(D, imageNo); % Update
				end
			end
		end
	end

	if options.makePlots==1
		if isfield(imgStats,'Area')==1
			openFigure(1996, '');
				subplot(2,1,1)
				hist(imgStats.Area,round(logspace(0,log10(max(imgStats.Area)))));
				box off;title('distribution of IC sizes');xlabel('area (px^2)');ylabel('count');
				set(gca,'xscale','log');
				h = findobj(gca,'Type','patch');
				set(h,'FaceColor',[0 0 0],'EdgeColor','w');
		end
		if ~isempty(options.valid)
			nPts = 2;
		else
			options.valid = ones(1,nImages);
			nPts = 1;
		end
		pointColors = ['g','r'];
		openFigure(1997, '');
			for pointNum = 1:nPts
				pointColor = pointColors(pointNum);
				if pointNum==1
					valid = logical(options.valid);
				else
					valid = logical(~options.valid);
				end
				openFigure(1997, '');
				fn=fieldnames(imgStats);
				for i=1:length(fn)
					subplot(2,ceil(length(fn)/2),i)
					% eval(['iStat=imgStats.' fn{i} ';']);
					% plot(find(valid),iStat(valid),[pointColor '.'])
					plot(find(valid),imgStats.(fn{i})(valid),[pointColor '.'])
					title(fn{i})
					hold on;box off;
					xlabel('rank'); ylabel(fn{i})
					hold off
				end
			end

	end
	function nUpdateParforProgress(~)
		if ~verLessThan('matlab', '9.2')
			p = p + 1;
			% if (mod(p,nInterval)==0||p==nImages)&&options_waitbarOn==1
			% 	cmdWaitbar(p,nImages,'','inputStr','','waitbarOn',1);
			% end
			if (mod(p,nInterval)==0||p==2||p==nImages)&&options_waitbarOn==1
				if p==nImages
					fprintf('%d\n',round(p/nImages*100))
				else
					fprintf('%d|',round(p/nImages*100))
				end
				% cmdWaitbar(p,nImages,'','inputStr','','waitbarOn',1);
			end
		end
	end
end

% OLD code

% regionStat
% figure;imagesc(inputImages(:,:,imageNo));
% for ifeature = featureList
% 	% regionStat = regionprops(iImage, ifeature{1});
% 	try
% 		% eval(['imgStats.' ifeature{1} '(imageNo) = regionStat.' ifeature{1} ';']);
% 		imgStats.(ifeature{1})(imageNo) = regionStat.(ifeature{1});
% 	catch
% 		% eval(['imgStats.' ifeature{1} '(imageNo) = NaN;']);
% 		imgStats.(ifeature{1})(imageNo) = NaN;
% 	end
% end

% regionStat = regionprops(iImage, 'Eccentricity','EquivDiameter','Area','Orientation','Perimeter','Solidity');
% imgStats.Eccentricity(imageNo) = regionStat.Eccentricity;
% imgStats.EquivDiameter(imageNo) = regionStat.EquivDiameter;
% imgStats.Area(imageNo) = regionStat.Area;
% imgStats.Orientation(imageNo) = regionStat.Orientation;
% imgStats.Perimeter(imageNo) = regionStat.Perimeter;
% imgStats.Solidity(imageNo) = regionStat.Solidity;

% if (imageNo==1||mod(imageNo,10)==0||imageNo==nImages)&&options.waitbarOn==1
% 	reverseStr = cmdWaitbar(imageNo,nImages,reverseStr,'inputStr','computing image features');
% end

% subplot(2,1,1)
% scatter3(imgStats.Eccentricity(valid),imgStats.Perimeter(valid),imgStats.Orientation(valid),[pointColor '.'])
% xlabel('Eccentricity');ylabel('perimeter');zlabel('Orientation');
% rotate3d on;hold on;
% subplot(2,1,2)
% scatter3(imgStats.Area(valid),imgStats.Perimeter(valid),imgStats.Solidity(valid),[pointColor '.'])
% xlabel('area');ylabel('perimeter');zlabel('solidity');
% rotate3d on;hold on;

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function jointfig(hfig,no_row,no_col)

%   jointfig(hfig,no_row,no_col)
%	--> joint subplots without any space between them
%   hfig : figure handler, if none, keyin gcf instead
%   no_row    : No. of row subplots
%   no_col    : No. of column subplots
%
%                 DO-SIG GONG,
%				  mail: D.Gong@soton.ac.uk
%				  HFRU/ISVR, University of Southampton, UK
%				  latest modified at 99-03-01 7:27PM

% All the movement of subplots should be done in unit of points
figure(hfig), hsubplot = get(hfig,'Children');
% convert the position unit from pixel into points : should be restored)
set(hfig,'unit','point')

% BEWARE! hsubplot has different order from the original subplot sequence
% for instance,
%
%  -----------------------         -----------------------
%  |     1    |     2     |        |     4    |     3     |
%  |----------+-----------|        |----------+-----------|
%  |     3    |     4     |        |     2    |     1     |
%  -----------------------         -----------------------
%       subplot(22i)                  get(gcf,'Children')
%
% THEREFORE, transpose hsubplot into the one in original subplot sequence, like this..

hsubplot = hsubplot(length(hsubplot):-1:1);

no_subplot1 = length(hsubplot);
no_space  = no_row*no_col;
no_delta = no_space - no_subplot1;

% in case of the odd number of subplots
if no_delta,
	for i = 1:no_delta
		addsubplot = subplot(no_row,no_col,no_subplot1+i);
		hsubplot = [hsubplot; addsubplot];
	end
end

no_subplot = length(hsubplot);

% Default position of figure in a window in point coord

for i=1:no_subplot,
	set(hsubplot(i),'unit','point'),
	tmp_ylab_pos = get(get(hsubplot(i),'ylabel'),'position');
	ylab_pos(i) = tmp_ylab_pos(1);
end

new_ylab_pos = min(ylab_pos);

coner1 = get(hsubplot(1),'position');
coner2 = get(hsubplot(length(hsubplot)),'position');

% position of lowest-left coner
inix = coner1(1);
iniy = coner2(2)*1.13;

% axis line width
alinewidth = get(hsubplot(1),'linewidth');

% total lengths
total_xlength = (coner2(1) + coner2(3) - coner1(1)) + (no_col-1) * alinewidth;
total_ylength = (coner1(2) + coner1(4) - coner2(2)) + (no_row-1) * alinewidth;

% width of each subplot
delx = 1.0 * total_xlength / no_col;

% height of each subplot
dely = 0.97 * total_ylength / no_row;

%index_loop = 0;                   % total subplots index
%for index_row = 1:no_row,         % loop for row index
%    for index_col = 1:no_col          % loop for column index
%        index_loop = index_loop+1;
%
%        startx = inix + (index_col - 1) * delx;
%        starty = iniy + (no_row - index_row) * dely;
%
%%.......It's kind of bug of MATLAB
%        if alinewidth < 1.0
%           set(hsubplot(index_loop),'position',...
%              [ startx - 0.5 * alinewidth * (index_col-1), ...
%                starty - 0.5 * alinewidth * (index_row-1), delx ,dely]);
%%              [startx-1.0*alinewidth*(index_col-1), starty+1.5*alinewidth*(index_row-1), delx ,dely]);
%        else
%           set(hsubplot(index_loop),'position',[startx, starty, delx ,dely]);
%        end
%
%        subplot(hsubplot(index_loop));

index_loop = no_subplot+1;              % total subplots index (reverse order)
for index_row = no_row:-1:1,             % loop for row index
    for index_col = no_col:-1:1          % loop for column index
        index_loop = index_loop - 1;

        startx = inix + (index_col - 1) * delx;
        starty = iniy + (no_row - index_row) * dely;
        POSITION = [startx, starty, delx ,dely];

%.......It's kind of bug of MATLAB
        if alinewidth < 1.0
           POSITION =  [ startx - 0.5 * alinewidth * (index_col-1), ...
                         starty + 0.9 * alinewidth * (index_row-1), delx ,dely];
%          POSITION =  [startx-1.0*alinewidth*(index_col-1), starty+1.5*alinewidth*(index_row-1), delx ,dely]);
        end

        set(hsubplot(index_loop),'position',POSITION);

        subplot(hsubplot(index_loop));

        iscale = size(get(gca,'yscale'),2);  % 3:log, 6:linear

        % remove xlabels & xticklabels of subplots located in upper rows other than lowest row

        if index_row ~= no_row,
        	if ~(no_delta & index_row == (no_row - 1) & index_col == no_col),
            	set(get(gca,'xlabel'),'String',[])
            	set(gca,'xticklabel',[]);  %remove xticklabel
            end
        end

        % remove ylabels & yticklabels of subplots located in right columns other than leftmost column
        if index_col ~= 1,
           set(get(gca,'ylabel'),'String',[])
           set(gca,'yticklabel',[]);  %remove yticklabel
        end

        % remove first yticklabel of subplots located in lower rows other than highest row, linear yscale only
		% .... only linear scale
        if index_row ~= 1 & iscale == 6
           	a = get(gca,'ytick'); b = get(gca,'ylim');
			if a(length(a)) == b(length(b)),
           		a = a(1:length(a)-1);
           		set(gca,'ytick',a);
           	end
        end

        % remove first xticklabel of subplots located in left columns other than rightmost column
		% .... only linear scale

		if ~no_delta,
  	      	if index_col ~= no_col & iscale == 6
       	       	a = get(gca,'xtick'); b = get(gca,'xlim');
           	   	if a(length(a)) == b(length(b)),
          			a = a(1:length(a)-1);
           			set(gca,'xtick',a);
           		end
        	end
        else
	   	    if index_col == no_col & index_row == no_row - 1 & iscale == 6,
       	       	a = get(gca,'xtick');
          		a = a(2:length(a));
           		set(gca,'xtick',a);
           	end
        end

   end
end

% get back to initial unit
set(hfig,'unit','default')
for i=1:no_subplot,	set(hsubplot(i),'unit','default'),end

% delete dummy subplots
if no_delta, for i = 1:no_delta, delete(hsubplot(no_subplot1+i)); end, end

author: Sable

repo: mcbench-benchmarks


text:
function [MStitch, result] = RGBMain_Process(MStitch, W_box, H_box, bdown)

[MStitch, result] = Fun_StitchRGB(MStitch.im2, W_box, H_box, bdown, MStitch);

author: decouples


text:
function [J, grad] = lrCostFunction(theta, X, y, lambda)
%LRCOSTFUNCTION Compute cost and gradient for logistic regression with
%regularization
%   J = LRCOSTFUNCTION(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Hint: The computation of the cost function and gradients can be
%       efficiently vectorized. For example, consider the computation
%
%           sigmoid(X * theta)
%
%       Each row of the resulting matrix will contain the value of the
%       prediction for that example. You can make use of this to vectorize
%       the cost function and gradient computations.
%
% Hint: When computing the gradient of the regularized cost function,
%       there're many possible vectorized solutions, but one solution
%       looks like:
%           grad = (unregularized gradient for logistic regression)
%           temp = theta;
%           temp(1) = 0;   % because we don't add anything for j = 0
%           grad = grad + YOUR_CODE_HERE (using the temp variable)
%

n = length(theta);

% X.shape = [m, n], y.shape = [m, 1], theta.shape = [n, 1]
% y_pred = sigmoid(X*theta), y_pred.shape = [m, 1]
y_pred = sigmoid(X*theta);
J = -1.0 / m * sum( y.*log(y_pred) + (1-y).*log(1-y_pred) );
J = J + lambda / (2.0*m) * sum( theta(2:n, 1).^2 );  % ignore theta(1)
grad = 1.0 / m .* (X' * (y_pred - y));
grad(2:n, 1) = grad(2:n, 1) + lambda / m * theta(2:n, 1);   % ignore theta(1)

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

grad = grad(:);

end

author: imLogM


text:
function [zipped,info] = norm2huff(vector)
%NORM2HUFF   Huffman codification (encoder)
%   For vectors, NORM2HUFF(X) returns a Huffman coded version of the input vector.
%   For matrices, X(:) is used as input.
%
%   Input must be of uint8 type, while the output is a uint8 array.
%
%   [...,INFO] = ... returns also a structure with data required to convert it back
%   to normal vector:
%
%      INFO.pad        = eventually added bits at the end of bit sequence;
%      INFO.huffcodes  = Huffman codewords;
%      INFO.ratio      = compression ratio;
%      INFO.length     = original data length;
%      INFO.maxcodelen = max codeword length;
%
%   Codewords are stored in the 52 available bits of a double. To avoid anbiguities,
%   after the last codeword bit, a "1" bit is added to terminate the codeword.
%   I.e. the max codeword length can be 51 bits.
%
%   See also HUFF2NORM

%   $Author: Giuseppe Ridino' $
%   $Revision: 1.0 $  $Date: 10-May-2004 15:03:04 $

% ensure to handle uint8 input vector
if ~isa(vector,'uint8'),
	error('input argument must be a uint8 vector')
end

% vector as a row
vector = vector(:)';

% frequency
f = frequency(vector);

% simbols presents in the vector are
simbols = find(f~=0); % first value is 1 not 0!!!
f = f(simbols);

% sort using the frequency
[f,sortindex] = sort(f);
simbols = simbols(sortindex);

% generate the codewords as the 52 bits of a double
len = length(simbols);
simbols_index = num2cell(1:len);
codeword_tmp = cell(len,1);
while length(f)>1,
	index1 = simbols_index{1};
	index2 = simbols_index{2};
	codeword_tmp(index1) = addnode(codeword_tmp(index1),uint8(0));
	codeword_tmp(index2) = addnode(codeword_tmp(index2),uint8(1));
	f = [sum(f(1:2)) f(3:end)];
	simbols_index = [{[index1 index2]} simbols_index(3:end)];
	% resort data in order to have the two nodes with lower frequency as first two
	[f,sortindex] = sort(f);
	simbols_index = simbols_index(sortindex);
end

% arrange cell array to have correspondance simbol <-> codeword
codeword = cell(256,1);
codeword(simbols) = codeword_tmp;

% calculate full string length
len = 0;
for index=1:length(vector),
	len = len+length(codeword{double(vector(index))+1});
end

% create the full 01 sequence
string = repmat(uint8(0),1,len);
pointer = 1;
for index=1:length(vector),
	code = codeword{double(vector(index))+1};
	len = length(code);
	string(pointer+(0:len-1)) = code;
	pointer = pointer+len;
end

% calculate if it is necessary to add padding zeros
len = length(string);
pad = 8-mod(len,8);
if pad>0,
	string = [string uint8(zeros(1,pad))];
end

% now save only usefull codewords
codeword = codeword(simbols);
codelen = zeros(size(codeword));
weights = 2.^(0:51);
maxcodelen = 0;
for index = 1:length(codeword),
	len = length(codeword{index});
	if len>maxcodelen,
		maxcodelen = len;
	end
	if len>0,
		code = sum(weights(codeword{index}==1));
		code = bitset(code,len+1);
		codeword{index} = code;
		codelen(index) = len;
	end
end
codeword = [codeword{:}];

% calculate zipped vector
cols = length(string)/8;
string = reshape(string,8,cols);
weights = 2.^(0:7);
zipped = uint8(weights*double(string));

% store data into a sparse matrix
huffcodes = sparse(1,1); % init sparse matrix
for index = 1:numel(codeword),
	huffcodes(codeword(index),1) = simbols(index);
end

% create info structure
info.pad = pad;
info.huffcodes = huffcodes;
info.ratio = cols./length(vector);
info.length = length(vector);
info.maxcodelen = maxcodelen;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function codeword_new = addnode(codeword_old,item)
codeword_new = cell(size(codeword_old));
for index = 1:length(codeword_old),
	codeword_new{index} = [item codeword_old{index}];
end


text:
function smallpot = marginalize_pot(bigpot, keep, maximize, useC)
% MARGINALIZE_POT Marginalize a mpot onto a smaller domain.
% smallpot = marginalize_pot(bigpot, keep, maximize, useC)
%
% The maximize argument is ignored - maxing out a Gaussian is the same as summing it out,
% since the mode and mean are equal.
% The useC argument is ignored.

node_sizes = sparse(1, max(bigpot.domain));
node_sizes(bigpot.domain) = bigpot.sizes;
sum_over = mysetdiff(bigpot.domain, keep);

[logp, mu, Sigma] = marginalize_gaussian(bigpot.logp, bigpot.mu, bigpot.Sigma, ...
					 keep, sum_over, node_sizes);
smallpot = mpot(keep, node_sizes(keep), logp, mu, Sigma);

%%%%%%

function [logpX, muX, SXX] = marginalize_gaussian(logp, mu, Sigma, X, Y, ns)
% MARGINALIZE_GAUSSIAN Compute Pr(X) from Pr(X,Y) where X and Y are jointly Gaussian.
% [logpX, muX, SXX] = marginalize_gaussian(logp, mu, Sigma, X, Y, ns)
%
% sizes(i) is the size of the i'th block in domain.

[muX, muY, SXX, SXY, SYX, SYY] = partition_matrix_vec(mu, Sigma, X, Y, ns);
logpX = logp; % Lauritzen (1996) p161

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function svm = svmRemove(svm, indices)
% svm = svmRemove(svm, indices)
% SUPPORT VECTOR MACHINE - REMOVE
% ---------------------------------------------------------
% Removes the specified condition from the SVM structure, patching up any
% holes it leaves in the order of the condition numbers.
%
% INPUTS
%   svm - Structure generated by svmInit.
%   indices - Vector of indices to remove from SVM structure.
%
% OUTPUTS
%   svm - Structure with specified indices removed.
%
% USAGE
%   After generating an SVM structure, removing conditions 2 and 3.
%       svm = svmInit(...);
%       svm = svmRemove(svm, [2 3]);
%
% See also SVMINIT, SVMRUN, SVMEXPORTMAP, SVMBOOTSTRAP, SVMRELABEL,
% SVMRUNOPTIONS, SLINIT.
%
% renobowen@gmail.com [2010]
%
    keepIndices = ~ismember(svm.group, indices);
    svm.run = svm.run(keepIndices, :);
    svm.group = svm.group(keepIndices, :);
    svm.trial = svm.trial(keepIndices, :);
    svm.data = svm.data(keepIndices, :);

    combineInto = min(svm.group);
    svm = svmRelabel(svm, [combineInto indices], svm.grouplabel{combineInto});
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% *************************************************************************************
% * DESCRIPTION:                                                                      *
% * This script computes the experiments performed in ref. [1].                       *
% * --------------------------------------------------------------------------------- *
% * REFERENCE:                                                                        *
% * [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and         *
% *     MRI texture features for the prediction of lung metastases in soft-tissue     *
% *     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),         *
% *     5471-5496. doi:10.1088/0031-9155/60/14/5471                                   *
% * --------------------------------------------------------------------------------- *
% * AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>                            *
% * --------------------------------------------------------------------------------- *
% * HISTORY:                                                                          *
% * - Creation: May 2015                                                              *
% * --------------------------------------------------------------------------------- *
% * STATEMENT:                                                                        *
% * This file is part of <https://github.com/mvallieres/radiomics/>,                  *
% * a package providing MATLAB programming tools for radiomics analysis.              *
% * --> Copyright (C) 2015  Martin Vallieres                                          *
% *                                                                                *
% *   This package is free software: you can redistribute it and/or modify            *
% *   it under the terms of the GNU General Public License as published by            *
% *   the Free Software Foundation, either version 3 of the License, or               *
% *   (at your option) any later version.                                             *
% *                                                                                *
% *   This package is distributed in the hope that it will be useful,                 *
% *   but WITHOUT ANY WARRANTY; without even the implied warranty of                  *
% *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                   *
% *   GNU General Public License for more details.                                    *
% *                                                                                *
% *   You should have received a copy of the GNU General Public License               *
% *   along with this package.  If not, see <http://www.gnu.org/licenses/>.           *
% *************************************************************************************

fprintf('\n')
help masterScript_STS
fprintf('\n')
warning off

% INITIALIZATION
pathWORK = pwd;
roiNumb = load('contour_Mass'); roiNumb = struct2cell(roiNumb); roiNumb = roiNumb{1};
outcome = load('outcome'); outcome = struct2cell(outcome); outcome = outcome{1};
nPatient = numel(outcome);

% TEXTURE EXTRACTION PARAMETERS AND DEGREES OF FREEDOM
MRIinv_cell = {'NoInv','Inv'};
MRIweight_mat = [1/4,1/3,1/2,2/3,3/4];
R_mat = [1/2,2/3,1,3/2,2];
scale_cell = {'pixelW',1,2,3,4,5};
algo_cell = {'Equal','Lloyd'};
Ng_mat = [8,16,32,64];
paramSEP = {R_mat,scale_cell,algo_cell,Ng_mat};                           % NOTE: paramSEP and paramFUS must always be of the same format and size for SEPARATE
paramFUS = {MRIinv_cell,MRIweight_mat,R_mat,scale_cell,algo_cell,Ng_mat}; % and FUSED SCANS as presented here, with the same ordering of different extraction parameters.
baselineSEP = [3,1,2,3]; % As defined in ref. [1]
baselineFUS = [1,3,3,1,2,3]; % As defined in ref. [1]
freedomSEP = [unique(perms([0 0 0 0]),'rows');unique(perms([0 0 0 1]),'rows');...
              unique(perms([0 0 1 1]),'rows');unique(perms([0 1 1 1]),'rows');...
              unique(perms([1 1 1 1]),'rows')];
freedomFUS = [unique(perms([0 0 0 0 0 0]),'rows');unique(perms([0 0 0 0 0 1]),'rows');...
              unique(perms([0 0 0 0 1 1]),'rows');unique(perms([0 0 0 1 1 1]),'rows');...
              unique(perms([0 0 1 1 1 1]),'rows');unique(perms([0 1 1 1 1 1]),'rows');...
              unique(perms([1 1 1 1 1 1]),'rows')];

% MULTIVARIABLE ANALYSIS PARAMETERS
nBoot = 1000;
alpha = 0.5; delta = 0.5;
setSize = 25;
fSetName = {'PET','SEPARATE','FUSED'};
maxOrder = 10;

% ************** START COMPUTATION **************

% 1. READ DATA DOWNLOADED FROM THE TCIA WEBSITE (http://dx.doi.org/10.7937/K9/TCIA.2015.7GO2GSKS)
fprintf('\n\n*********************** ORGANIZING AND PROCESSING DICOM DATA FROM TCIA WEBSITE ***********************\n')
readAllDICOM_STS([pathWORK,'/Soft-tissue-Sarcoma'],nPatient)

% 2. COMPUTE NON-TEXTURE FEATURES
fprintf('\n\n*********************** COMPUTING NON-TEXTURE FEATURES ***********************\n')
calcAllNonTextureFeatures_STS(pathWORK,nPatient,roiNumb,outcome)

% 3. COMPUTE AND ORGANIZING ALL TEXTURE FEATURES
mkdir('TEXTURES'), fprintf('\n')
calcAllSeparateTextures_STS(pathWORK,nPatient,roiNumb,R_mat,scale_cell,algo_cell,Ng_mat)
organizeSeparateTextures_STS(pathWORK,nPatient,R_mat,scale_cell,algo_cell,Ng_mat)
calcAllFusedTextures_STS(pathWORK,nPatient,roiNumb,MRIinv_cell,MRIweight_mat,R_mat,scale_cell,algo_cell,Ng_mat)
organizeFusedTextures_STS(pathWORK,nPatient,MRIinv_cell,MRIweight_mat,R_mat,scale_cell,algo_cell,Ng_mat)

% 5. PERFORM FEATURE SET REDUCTION
mkdir('FSET')
nonTextures = load('nonTextures'); nonTextures = struct2cell(nonTextures); nonTextures = nonTextures{1};
PET  = load('textures_PET');  PET  = struct2cell(PET);  PET  = PET{1};
T1   = load('textures_T1');   T1   = struct2cell(T1);   T1   = T1{1};
T2FS = load('textures_T2FS'); T2FS = struct2cell(T2FS); T2FS = T2FS{1};
PET_T1   = load('textures_PET_T1');   PET_T1   = struct2cell(PET_T1);   PET_T1   = PET_T1{1};
PET_T2FS = load('textures_PET_T2FS'); PET_T2FS = struct2cell(PET_T2FS); PET_T2FS = PET_T2FS{1};

tic
fprintf('\n\nFinding the path to the ''MINE.jar'' application on the system ... ')
pathMINE = findMINE('Linux');
fprintf('DONE\n')
toc

% For Feature set 1: PET
fprintf('\n*********************** PERFORMING FEATURE SET REDUCTION FOR ''PET'' FEATURE SET ***********************\n')
calcAllFeatureSets_STS(pathWORK,pathMINE,fSetName{1},outcome,setSize,nonTextures,{PET},{'PET'},paramSEP,freedomSEP,baselineSEP,alpha,delta,nBoot)

% For Feature set 2: PET, T1, T2FS
fprintf('\n*********************** PERFORMING FEATURE SET REDUCTION FOR ''SEPARATE'' FEATURE SET **********************\n')
calcAllFeatureSets_STS(pathWORK,pathMINE,fSetName{2},outcome,setSize,nonTextures,{PET,T1,T2FS},{'PET','T1','T2FS'},paramSEP,freedomSEP,baselineSEP,alpha,delta,nBoot)

% For Feature set 3: PET_T1, PET_T2FS
fprintf('\n*********************** PERFORMING FEATURE SET REDUCTION FOR ''FUSED'' FEATURE SET ***********************\n')
calcAllFeatureSets_STS(pathWORK,pathMINE,fSetName{3},outcome,setSize,nonTextures,{PET_T1,PET_T2FS},{'PET_T1','PET_T2FS'},paramFUS,freedomFUS,baselineFUS,alpha,delta,nBoot)

% 6. PERFORM FEATURE SET SELECTION
mkdir('MODELS')
fprintf('\n\nFinding the path to''fastAUC.cpp'' on the system --> COMPILATION ... ')
try
    compileFastAUC('Linux')
    fprintf('DONE\n')
catch
    fprintf('FAILED (AUC computations will be slower)\n')
end

% For Feature set 1: PET
fprintf('\n*********************** PERFORMING FEATURE SELECTION FOR ''PET'' FEATURE SET ***********************\n')
computeAllModelChoice_STS(pathWORK,fSetName{1},outcome,freedomSEP,maxOrder,nBoot)

% For Feature set 2: PET, T1, T2FS
fprintf('\n*********************** PERFORMING FEATURE SELECTION  FOR ''SEPARATE'' FEATURE SET ***********************\n')
computeAllModelChoice_STS(pathWORK,fSetName{2},outcome,freedomSEP,maxOrder,nBoot)

% For Feature set 3: PET_T1, PET_T2FS
fprintf('\n*********************** PERFORMING FEATURE SELECTION  FOR ''FUSED'' FEATURE SET ***********************\n')
computeAllModelChoice_STS(pathWORK,fSetName{3},outcome,freedomFUS,maxOrder,nBoot)

% 7. PERFORM PREDICTION PERFORMANCE ESTIMATION
mkdir('RESULTS'), fprintf('\n')

% For Feature set 1: PET
fprintf('\n*********************** PERFORMING PREDICTION PERFORMANCE ESTIMATION FOR ''PET'' FEATURE SET ***********************\n')
computeAllPrediction_STS(pathWORK,fSetName{1},outcome,freedomSEP,maxOrder,nBoot)

% For Feature set 2: PET, T1, T2FS
fprintf('\n*********************** PERFORMING PREDICTION PERFORMANCE ESTIMATION FOR ''SEPARATE'' FEATURE SET ***********************\n')
computeAllPrediction_STS(pathWORK,fSetName{2},outcome,freedomSEP,maxOrder,nBoot)

% For Feature set 3: PET_T1, PET_T2FS
fprintf('\n*********************** PERFORMING PREDICTION PERFORMANCE ESTIMATION FOR ''FUSED'' FEATURE SET ***********************\n')
computeAllPrediction_STS(pathWORK,fSetName{3},outcome,freedomFUS,maxOrder,nBoot)

% Finding the best combinations of model order and texture extraction parameter degree of freedom for all feature set types
groupExperiments_STS(pathWORK,fSetName{1},freedomSEP,maxOrder)
groupExperiments_STS(pathWORK,fSetName{2},freedomSEP,maxOrder)
groupExperiments_STS(pathWORK,fSetName{3},freedomFUS,maxOrder)

% 8. CHOICE OF BEST PARSIMONIOUS MODEL (requires user input)
fprintf('\n')
plotPredictionResults_STS([pathWORK,'/RESULTS'],fSetName,{'AUC632','Sensitivity632','Specificity632'},maxOrder)
while 1
    set = input(['\nWhich feature set provides the best parsimonious model? \n' ...
                 '--> For the ',fSetName{1},' feature set, type ''1'' and press ENTER \n' ...
                 '--> For the ',fSetName{2},' feature set, type ''2'' and press ENTER \n' ...
                 '--> For the ',fSetName{3},' feature set, type ''3'' and press ENTER \n' ...
                 'ANSWER: ']);
    fprintf('\n')
    if isnumeric(set) && (set == 1 || set == 2 || set == 3)
        break
    end
end
while 1
    order = input(['Which model order of the ',fSetName{set},' feature set provides the best parsimonious model? \n' ...
                   '--> Type a number between 1 to ',num2str(maxOrder),' and press ENTER \n' ...
                   'ANSWER: ']);
    fprintf('\n')
    if isnumeric(order) && order <= 10 && order >= 1
        break
    end
end
cd([pathWORK,'/RESULTS'])
results = load(['RESULTS_',fSetName{set},'_BEST']); results = struct2cell(results); results = results{1};
finalModel = results.(['Order',num2str(order)]);
cd(pathWORK), mkdir('FINAL_MODEL'), cd('FINAL_MODEL'), save('finalModel', 'finalModel')

% 9. COMPUTING THE LOGISTIC REGRESSION COEFFICIENTS AND BOOTSTRAP CONFIDENCE INTERVALS OF THE FINAL MODEL
fprintf('\n\nCOMPUTING THE LOGISTIC REGRESSION COEFFICIENTS OF THE FINAL MODEL ... ')
tic
[coeff,response,modelCI] = computeModelCoefficients(finalModel.Data,outcome,'IABR');
save('coeff','coeff'), save('response','response'), save('modelCI','modelCI')
fprintf('DONE\n')
toc

% 10. DISPLAYING THE FINAL MODEL AND CORRESPONDING PREDICTION PERFORMANCE ESTIMATION
plotSigmoidalResponse(response,outcome,modelCI,'LungMets')
fprintf(['\n\n\n --> THE FINAL MULTIVARIABLE MODEL IS:\n\n'...
         '               g(x) =               \n'])
for i = 1:order
    fprintf([num2str(coeff(i)),' X ',finalModel.Name{i},'\n'])
    fprintf('                    +               \n')
end
fprintf(['                   ',num2str(coeff(end)),'\n'])
fprintf('\nWITH CORRESPONDING PREDICTION PERFORMANCE ESTIMATION:\n')
fprintf(['AUC = ',num2str(roundsd(finalModel.AUC632,ceil(log10(finalModel.AUC632/roundsd(finalModel.SE_AUC632,1))))),' ± ',num2str(roundsd(finalModel.SE_AUC632,1)),'\n'])
fprintf(['Sensitivity = ',num2str(roundsd(finalModel.Sensitivity632,ceil(log10(finalModel.Sensitivity632/roundsd(finalModel.SE_Sensitivity632,1))))),' ± ',num2str(roundsd(finalModel.SE_Sensitivity632,1)),'\n'])
fprintf(['Specificity = ',num2str(roundsd(finalModel.Specificity632,ceil(log10(finalModel.Specificity632/roundsd(finalModel.SE_Specificity632,1))))),' ± ',num2str(roundsd(finalModel.SE_Specificity632,1)),'\n'])
fprintf('\n')
cd(pathWORK)

author: mvallieres

repo: radiomics


text:
function setG(OCP,G)
    % assert
    if OCP.dim.u + OCP.dim.x < 30
        Gux = jacobian(G,[OCP.u;OCP.x]);
        hasUX = has(Gux(:),[OCP.u;OCP.x]);
        if sum(hasUX(:)) ~= 0
            error('G must be a linear function of u and x!');
        end
    end
    OCP.G = G;
    % Global variable
    global ParNMPCGlobalVariable
    G_formula = formula(OCP.G);
    [zDim,~] = size(G_formula);
    ParNMPCGlobalVariable.dim.z = zDim;
    ParNMPCGlobalVariable.solutionInitGuess.z = zeros(zDim,ParNMPCGlobalVariable.N);
    OCP.dim.z      = zDim;
    OCP.z          = sym('z',[OCP.dim.z,1]);
    if size(OCP.z,1) ~= OCP.dim.z
       OCP.z = OCP.z.';
    end
    OCP.LBarrier   = sym(0);
    %% barrier term
    for i = 1:zDim
        OCP.LBarrier = OCP.LBarrier - log(G_formula(i));
    end
    %% linear damping term
    for i = 1:zDim
        OCP.LBarrier = OCP.LBarrier + 1e-4*G_formula(i);
    end
end

author: deng-haoyang

repo: ParNMPC


text:
% Analytic SME meanamp calc snippet
% axs Dec 2018
% ams Jan 2023, includes corrected SME
function [outdata] = sme_analytic(EEG, epoch_list, dq_window_times)

sizes = size(EEG.data);

assert(numel(sizes)==3, 'sme_analytic needs bin-epoched EEG data');

if exist('dq_window_times','var') == 0
    dq_window_times = [];
end

% if there are 3 cols to win_times, use last 2
if size(dq_window_times,2) == 3
    dq_window_times = dq_window_times(:,2:3);
end

n_bins = length(epoch_list);
n_elec = sizes(1);
n_times = sizes(2);
n_beps_total = sizes(3);

if isempty(dq_window_times)
    n_windows = 1;
    win_times_starts = 1;
    win_times_ends = EEG.times(end);
else
    n_windows = size(dq_window_times,1);
    win_times_starts = dq_window_times(:,1);
    win_times_ends = dq_window_times(:,2);
end

% Convert times from ms to datapoint idx
win_dps_starts = zeros(1,n_windows);
win_dps_ends = zeros(1,n_windows);
for t = 1:n_windows
    win_dps_starts(t) = find(abs(EEG.times-win_times_starts(t))==min(abs(EEG.times-win_times_starts(t))),1);
    win_dps_ends(t) = find(abs(EEG.times-win_times_ends(t))==min(abs(EEG.times-win_times_ends(t))),1);
end

for b = 1:n_bins
    n_beps_per_bin(b) = numel(epoch_list(b).good_bep_indx);
end

SME_out = zeros(n_elec,n_windows,n_bins);
SME_unbias_out = zeros(n_elec,n_windows,n_bins);

for b = 1:n_bins        % for each bin
    for t = 1:n_windows % and each SME time window

    data_here = EEG.data(1:n_elec,win_dps_starts(t):win_dps_ends(t),epoch_list(b).good_bep_indx);
    window_mean = squeeze(mean(data_here,2));

    %SD (N-1)
    SME_sd = std(window_mean,0,2);
    %SME_bias = std(window_mean,0,2) / sqrt(n_beps_per_bin(b));
    SME_bias = SME_sd / sqrt(n_beps_per_bin(b));
    SME_out(1:n_elec,t,b) = SME_bias;

    %Unbiased SD (Gurland & Tripathi)
    nTimes = size(window_mean,2);
    SME_unbias = (((SME_sd * sqrt((nTimes-1)))/ sqrt((nTimes-(3/2)+(1/(8*(nTimes-1)))))) / sqrt(n_beps_per_bin(b)));
    SME_unbias_out(1:n_elec,t,b) = SME_unbias;

    end
end
%combine SME structs
outdata = struct();
outdata.SME = SME_out;
outdata.SME_corr = SME_unbias_out;

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
%%*****************************************************************
%% NTpred: Compute (dX,dy,dZ) for NT direction.
%%
%% compute SVD of Xchol*Zchol via eigenvalue decompostion of
%%     Zchol * X * Zchol' = V * diag(sv2) * V'.
%% compute W satisfying W*Z*W = X.
%%     W = G'*G,  where G = diag(sqrt(sv)) * (invZchol*V)'
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

 function [par,dX,dy,dZ,coeff,L,hRd] = ...
          NTpred(blk,At,par,rp,Rd,sigmu,X,Z,Zchol,invZchol);

    global schurfun schurfun_par
%%
%% compute NT scaling matrix
%%
    [par.W,par.G,par.sv,par.gamx,par.gamz,par.dd,par.ee,par.ff] = ...
     NTscaling(blk,X,Z,Zchol,invZchol);
%%
%% compute schur matrix
%%
    m = length(rp);
    schur = sparse(m,m);
    UU = []; EE = []; Afree = [];
    dX = cell(size(blk,1),1); dy = []; dZ = cell(size(blk,1),1);
%%
    for p = 1:size(blk,1)
       pblk = blk(p,:);
       if strcmp(pblk{1},'l')
          [schur,UU,EE] = schurmat_lblk(blk,At,par,schur,UU,EE,p,par.dd);
       elseif strcmp(pblk{1},'q');
          [schur,UU,EE] = schurmat_qblk(blk,At,par,schur,UU,EE,p,par.dd,par.ee);
       elseif strcmp(pblk{1},'s')
          if isempty(schurfun{p})
             schur = schurmat_sblk(blk,At,par,schur,p,par.W);
          elseif isstr(schurfun{p})
             schurtmp = sparse(m,m);
             if ~isempty(par.permZ{p})
                Wp = par.W{p}(par.permZ{p},par.permZ{p});
             else
                Wp = par.W{p};
             end
             eval(['schurtmp = ',schurfun{p},'(Wp,Wp,schurfun_par(p,:));']);
             schur = schur + schurtmp;
          end
       elseif strcmp(pblk{1},'u')
          Afree = [Afree, At{p}'];
       end
    end
%%
%% compute rhs
%%
    [rhs,EinvRc,hRd] = NTrhsfun(blk,At,par,X,Z,rp,Rd,sigmu);
%%
%% solve linear system
%%
    [xx,coeff,L] = linsysolve(par,schur,UU,Afree,EE,rhs);
%%
%% compute (dX,dZ)
%%
    [dX,dy,dZ] = NTdirfun(blk,At,par,Rd,EinvRc,xx,m);
%%**********************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function varargout = subsref(varargin)
%SUBSREF   CHEBFUN2 subsref.
% ( )
%   F(X, Y) returns the values of the CHEBFUN2 F evaluated at (X,Y). See
%   CHEBFUN/FEVAL for further details. F(:, Y) returns a chebfun
%   representing the function F along that column slice, and F(X, :)
%   returns a chebfun representing F along that row slice. F(:, :) returns
%   F.
%
%   F(G) computes the composition of F and G where G is a CHEBFUN with two
%   columns, a CHEBFUN2V, a CHEBFUN3V with two components, or a DISKFUNV.
%   If G is a CHEBFUN with one column, a CHEBFUN2, a CHEBFUN3, a DISKFUN or
%   a SPHEREFUN, F(G) is interpreted as F(real(G), imag(G)), regardless of
%   whether G is real or complex.
%
%   F(X, Y) with CHEBFUNs X and Y returns the CHEBFUN G(t) = F(X(t), Y(t)).
%   If X and Y are CHEBFUN2 objects, then F(X, Y) is a CHEBFUN2.
%   If X and Y are CHEBFUN3 objects, then F(X, Y) is a CHEBFUN3.
%
% .
%   F.PROP returns the property PROP of F as defined by GET(F, 'PROP').
%
% {}
%   F{S1, S2, S3, S4} restricts F to the domain [S1, S2, S3, S4]. See
%   CHEBFUN2/RESTRICT for further details. Note that F{[S1,S2, S3, S4]} is
%   not supported due to the behaviour of the MATLAB subsref() command.
%
% See also FEVAL, GET, RESTRICT, SUBSREF.

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

[varargout{1:nargout}] = subsref@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [contact, fric_coef, geometry] = RightSoleInside(robot)

    param = sys.GetExtraParams();

    r_foot_frame = robot.Joints(getJointIndices(robot, 'r_leg_akx'));
    contact = CoordinateFrame(...
        'Name','RightSole',...
        'Reference',r_foot_frame,...
        'Offset',[0, param.wf/2, param.hf],...
        'R',[0,0,0]... % z-axis is the normal axis, so no rotation required
        );

    fric_coef.mu = param.mu;
    fric_coef.gamma = param.gamma;

    geometry.la = param.wf/2;
    geometry.lb = param.wf/2;
    geometry.La = param.lt;
    geometry.Lb = param.lh;
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
classdef Step1 < adi.StepBase

    %   Copyright 2017 The MathWorks, Inc.
    properties (Hidden)
       InstallButton
    end

    methods
        function obj = Step1(varargin)
            obj@adi.StepBase(varargin{:});

            %       obj.Workflow.InstallDriverWorkflow = true;

            parentContainer = obj.ContentPanel;

            obj.Title.Text = 'ADI Board Support Package Installer';

            obj.InfoText1.Position = [20 340 400 30];
            obj.InfoText1.Text = 'Step 1: Installing drivers';

            obj.InfoText2.Position = [20 290 400 50];
            obj.InfoText2.Text = 'Test2';

            obj.InstallButton = matlab.hwmgr.internal.hwsetup.Button.getInstance(parentContainer);
            obj.InstallButton.ButtonPushedFcn = @obj.LinuxInstall;
            obj.InstallButton.Text = 'WINNER';

            obj.HelpText.AboutSelection = 'About Selection';
            obj.HelpText.WhatToConsider = 'Test3';
        end

        function out = getNextScreenID(obj)
            % HSA infrastructure requires a "getNextScreenID" method implemented
            % in the leaf class. Otherwise, it renders "Finish" button instead of
            % "Next".
            out = getNextScreenIDImpl(obj);
        end

        function restoreScreen(obj)
            obj.enableScreen();
        end
    end

    methods (Access = protected)

        function out = getNextScreenIDImpl(~)
            out = 'adi.Step2';
        end

        function LinuxInstall(~,~,~)
            system('sudo whoami');
        end

    end

end

author: analogdevicesinc


text:
%isimatrix - Test if parameter is a matrix of integers (>= 2 columns).
%
%  USAGE
%
%    test = isimatrix(x,test1,test2,...)
%
%    x              parameter to test
%    test1...       optional list of additional tests
%
%  EXAMPLES
%
%    % Test if x is a matrix of doubles
%    isimatrix(x)
%
%    % Test if x is a matrix of strictly positive doubles
%    isimatrix(x,'>0')
%
%  NOTE
%
%    The tests ignore NaNs, e.g. isimatrix([500 nan;4 79]), isimatrix([1 nan 3],'>0') and
%    isimatrix([nan -7;nan nan;-2 -5],'<=0') all return 1.
%
%  SEE ALSO
%
%    See also isdmatrix, isdvector, isdscalar, isivector, isiscalar, isstring,
%    islscalar, islvector, islmatrix.
%

% Copyright (C) 2010 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

function test = isimatrix(x,varargin)

% Check number of parameters
if nargin < 1,
  error('Incorrect number of parameters (type ''help isimatrix'' for details).');
end

% Test: doubles, two dimensions, two or more columns?
test = isa(x,'double') & length(size(x)) == 2 & size(x,2) >= 2;

% Ignore NaNs (this reshapes the matrix, but it does not matter for the remaining tests)
x = x(~isnan(x));

% Test: integers?
test = test & all(round(x)==x);

% Optional tests
for i = 1:length(varargin),
	try
		if ~eval(['all(x(:)' varargin{i} ');']), test = false; return; end
	catch err
		error(['Incorrect test ''' varargin{i} ''' (type ''help isimatrix'' for details).']);
	end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip