text:
function [SPM,xSPM] = spm_getSPM(varargin)
% Compute a specified and thresholded SPM/PPM following estimation
% FORMAT [SPM,xSPM] = spm_getSPM;
% Query SPM in interactive mode.
%
% FORMAT [SPM,xSPM] = spm_getSPM(xSPM);
% Query SPM in batch mode. See below for a description of fields that may
% be present in xSPM input. Values for missing fields will be queried
% interactively.
%
% xSPM      - structure containing SPM, distribution & filtering details
% .swd      - SPM working directory - directory containing current SPM.mat
% .title    - title for comparison (string)
% .Z        - minimum of Statistics {filtered on u and k}
% .n        - conjunction number <= number of contrasts
% .STAT     - distribution {Z, T, X, F or P}
% .df       - degrees of freedom [df{interest}, df{residual}]
% .STATstr  - description string
% .Ic       - indices of contrasts (in SPM.xCon)
% .Im       - indices of masking contrasts (in xCon)
% .pm       - p-value for masking (uncorrected)
% .Ex       - flag for exclusive or inclusive masking
% .u        - height threshold
% .k        - extent threshold {voxels}
% .XYZ      - location of voxels {voxel coords}
% .XYZmm    - location of voxels {mm}
% .S        - search Volume {voxels}
% .R        - search Volume {resels}
% .FWHM     - smoothness {voxels}
% .M        - voxels -> mm matrix
% .iM       - mm -> voxels matrix
% .VOX      - voxel dimensions {mm} - column vector
% .DIM      - image dimensions {voxels} - column vector
% .Vspm     - Mapped statistic image(s)
% .Ps       - uncorrected P values in searched volume (for voxel FDR)
% .Pp       - uncorrected P values of peaks (for peak FDR)
% .Pc       - uncorrected P values of cluster extents (for cluster FDR)
% .uc       - 0.05 critical thresholds for FWEp, FDRp, FWEc, FDRc
% .thresDesc - description of height threshold (string)
%
% Required fields of SPM
%
% xVol   - structure containing details of volume analysed
%
% xX     - Design Matrix structure
%        - (see spm_spm.m for structure)
%
% xCon   - Contrast definitions structure array
%        - (see also spm_FcUtil.m for structure, rules & handling)
% .name  - Contrast name
% .STAT  - Statistic indicator character ('T', 'F' or 'P')
% .c     - Contrast weights (column vector contrasts)
% .X0    - Reduced design matrix data (spans design space under Ho)
%          Stored as coordinates in the orthogonal basis of xX.X from spm_sp
%          (Matrix in SPM99b)  Extract using X0 = spm_FcUtil('X0',...
% .iX0   - Indicates how contrast was specified:
%          If by columns for reduced design matrix then iX0 contains the
%          column indices. Otherwise, it's a string containing the
%          spm_FcUtil 'Set' action: Usually one of {'c','c+','X0'}
% .X1o   - Remaining design space data (X1o is orthogonal to X0)
%          Stored as coordinates in the orthogonal basis of xX.X from spm_sp
%          (Matrix in SPM99b)  Extract using X1o = spm_FcUtil('X1o',...
% .eidf  - Effective interest degrees of freedom (numerator df)
%        - Or effect-size threshold for Posterior probability
% .Vcon  - Name of contrast (for 'T's) or ESS (for 'F's) image
% .Vspm  - Name of SPM image
%
% Evaluated fields in xSPM (input)
%
% xSPM      - structure containing SPM, distribution & filtering details
% .swd      - SPM working directory - directory containing current SPM.mat
% .title    - title for comparison (string)
% .Ic       - indices of contrasts (in SPM.xCon)
% .n        - conjunction number <= number of contrasts
% .Im       - indices of masking contrasts (in xCon)
% .pm       - p-value for masking (uncorrected)
% .Ex       - flag for exclusive or inclusive masking
% .u        - height threshold
% .k        - extent threshold {voxels}
% .thresDesc - description of height threshold (string)
%
% In addition, the xCon structure is updated. For newly evaluated
% contrasts, SPM images (spmT_????.{img,hdr}) are written, along with
% contrast (con_????.{img,hdr}) images for SPM{T}'s, or Extra
% Sum-of-Squares images (ess_????.{img,hdr}) for SPM{F}'s.
%
% The contrast images are the weighted sum of the parameter images,
% where the weights are the contrast weights, and are uniquely
% estimable since contrasts are checked for estimability by the
% contrast manager. These contrast images (for appropriate contrasts)
% are suitable summary images of an effect at this level, and can be
% used as input at a higher level when effecting a random effects
% analysis. (Note that the ess_????.{img,hdr} and
% SPM{T,F}_????.{img,hdr} images are not suitable input for a higher
% level analysis.) See spm_RandFX.man for further details.
%
%__________________________________________________________________________
%
% spm_getSPM prompts for an SPM and applies thresholds {u & k}
% to a point list of voxel values (specified with their locations {XYZ})
% This allows the SPM be displayed and characterized in terms of regionally
% significant effects by subsequent routines.
%
% For general linear model Y = XB + E with data Y, design matrix X,
% parameter vector B, and (independent) errors E, a contrast c'B of the
% parameters (with contrast weights c) is estimated by c'b, where b are
% the parameter estimates given by b=pinv(X)*Y.
%
% Either single contrasts can be examined or conjunctions of different
% contrasts. Contrasts are estimable linear combinations of the
% parameters, and are specified using the SPM contrast manager
% interface [spm_conman.m]. SPMs are generated for the null hypotheses
% that the contrast is zero (or zero vector in the case of
% F-contrasts). See the help for the contrast manager [spm_conman.m]
% for a further details on contrasts and contrast specification.
%
% A conjunction assesses the conjoint expression of multiple effects. The
% conjunction SPM is the minimum of the component SPMs defined by the
% multiple contrasts.  Inference on the minimum statistics can be
% performed in different ways.  Inference on the Conjunction Null (one or
% more of the effects null) is accomplished by assessing the minimum as
% if it were a single statistic; one rejects the conjunction null in
% favor of the alternative that k=nc, that the number of active effects k
% is equal to the number of contrasts nc.  No assumptions are needed on
% the dependence between the tests.
%
% Another approach is to make inference on the Global Null (all effects
% null).  Rejecting the Global Null of no (u=0) effects real implies an
% alternative that k>0, that one or more effects are real.   A third
% Intermediate approach, is to use a null hypothesis of no more than u
% effects are real.  Rejecting the intermediate null that k<=u implies an
% alternative that k>u, that more than u of the effects are real.
%
% The Global and Intermediate nulls use results for minimum fields which
% require the SPMs to be identically distributed and independent. Thus,
% all component SPMs must be either SPM{t}'s, or SPM{F}'s with the same
% degrees of freedom. Independence is roughly guaranteed for large
% degrees of freedom (and independent data) by ensuring that the
% contrasts are "orthogonal". Note that it is *not* the contrast weight
% vectors per se that are required to be orthogonal, but the subspaces of
% the data space implied by the null hypotheses defined by the contrasts
% (c'pinv(X)). Furthermore, this assumes that the errors are
% i.i.d. (i.e. the estimates are maximum likelihood or Gauss-Markov. This
% is the default in spm_spm).
%
% To ensure approximate independence of the component SPMs in the case of
% the global or intermediate null, non-orthogonal contrasts are serially
% orthogonalised in the order specified, possibly generating new
% contrasts, such that the second is orthogonal to the first, the third
% to the first two, and so on.  Note that significant inference on the
% global null only allows one to conclude that one or more of the effects
% are real.  Significant inference on the conjunction null allows one to
% conclude that all of the effects are real.
%
% Masking simply eliminates voxels from the current contrast if they
% do not survive an uncorrected p value (based on height) in one or
% more further contrasts.  No account is taken of this masking in the
% statistical inference pertaining to the masked contrast.
%
% The SPM is subject to thresholding on the basis of height (u) and the
% number of voxels comprising its clusters {k}. The height threshold is
% specified as above in terms of an [un]corrected p value or
% statistic.  Clusters can also be thresholded on the basis of their
% spatial extent. If you want to see all voxels simply enter 0.  In this
% instance the 'set-level' inference can be considered an 'omnibus test'
% based on the number of clusters that obtain.
%
% BAYESIAN INFERENCE AND PPMS - POSTERIOR PROBABILITY MAPS
%
% If conditional estimates are available (and your contrast is a T
% contrast) then you are asked whether the inference should be 'Bayesian'
% or 'classical' (using GRF).  If you choose Bayesian the contrasts are of
% conditional (i.e. MAP) estimators and the inference image is a
% posterior probability map (PPM).  PPMs encode the probability that the
% contrast exceeds a specified threshold.  This threshold is stored in
% the xCon.eidf.  Subsequent plotting and tables will use the conditional
% estimates and associated posterior or conditional probabilities.
%
% see spm_results_ui.m for further details of the SPM results section.
% see also spm_contrasts.m
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes, Karl Friston & Jean-Baptiste Poline
% $Id$

%-GUI setup
%--------------------------------------------------------------------------
spm_help('!ContextHelp',mfilename)
spm('Pointer','Arrow')

%-Select SPM.mat & note SPM results directory
%--------------------------------------------------------------------------
if (nargin > 0)
    xSPM = varargin{1};
end
try
    swd = xSPM.swd;
    sts = 1;
catch
    [spmmatfile, sts] = spm_select(1,'^SPM\.mat$','Select SPM.mat');
    swd = spm_str_manip(spmmatfile,'H');
end;
if ~sts, SPM = []; xSPM = []; return; end

%-Preliminaries...
%==========================================================================

%-Load SPM.mat
%--------------------------------------------------------------------------
try
    load(fullfile(swd,'SPM.mat'));
catch
    error(['Cannot read ' fullfile(swd,'SPM.mat')]);
end
SPM.swd = swd;

%-Change directory so that relative filenames are valid
%--------------------------------------------------------------------------
cd(SPM.swd);

%-Check the model has been estimated
%--------------------------------------------------------------------------
try
    XYZ  = SPM.xVol.XYZ;
catch

    %-Check the model has been estimated
    %----------------------------------------------------------------------
    str = { 'This model has not been estimated.';...
        'Would you like to estimate it now?'};
    if spm_input(str,1,'bd','yes|no',[1,0],1)
        SPM = spm_spm(SPM);
    else
        SPM = []; xSPM = [];
        return
    end
end

xX   = SPM.xX;                      %-Design definition structure
XYZ  = SPM.xVol.XYZ;                %-XYZ coordinates
S    = SPM.xVol.S;                  %-search Volume {voxels}
R    = SPM.xVol.R;                  %-search Volume {resels}
M    = SPM.xVol.M(1:3,1:3);         %-voxels to mm matrix
VOX  = sqrt(diag(M'*M))';           %-voxel dimensions

% check the data and other files have valid filenames
%-----------------------------------------------------------------------
try, SPM.xY.VY     = spm_check_filename(SPM.xY.VY);     end
try, SPM.xVol.VRpv = spm_check_filename(SPM.xVol.VRpv); end
try, SPM.Vbeta     = spm_check_filename(SPM.Vbeta);     end
try, SPM.VResMS    = spm_check_filename(SPM.VResMS);    end
try, SPM.VM        = spm_check_filename(SPM.VM);        end

%-Contrast definitions
%==========================================================================

%-Load contrast definitions (if available)
%--------------------------------------------------------------------------
try
    xCon = SPM.xCon;
catch
    xCon = {};
end

%==========================================================================
% - C O N T R A S T S ,   S P M    C O M P U T A T I O N ,    M A S K I N G
%==========================================================================

%-Get contrasts
%--------------------------------------------------------------------------
try
    Ic = xSPM.Ic;
catch
    [Ic,xCon] = spm_conman(SPM,'T&F',Inf,...
                           '    Select contrasts...',' for conjunction',1);
end
if isempty(xCon)
    % figure out whether new contrasts were defined, but not selected
    % do this by comparing length of SPM.xCon to xCon, remember added
    % indices to run spm_contrasts on them as well
    try
        noxCon = numel(SPM.xCon);
    catch
        noxCon = 0;
    end
    IcAdd = (noxCon+1):numel(xCon);
else
    IcAdd = [];
end

nc        = length(Ic);  % Number of contrasts

%-Allow user to extend the null hypothesis for conjunctions
%
% n: conjunction number
% u: Null hyp is k<=u effects real; Alt hyp is k>u effects real
%    (NB Here u is from Friston et al 2004 paper, not statistic thresh).
%                  u         n
% Conjunction Null nc-1      1     |    u = nc-n
% Intermediate     1..nc-2   nc-u  |    #effects under null <= u
% Global Null      0         nc    |    #effects under alt  > u,  >= u+1
%----------------------------------+---------------------------------------
if (nc > 1)
    try
        n = xSPM.n;
    catch
        if nc==2
            But = 'Conjunction|Global';      Val=[1 nc];
        else
            But = 'Conj''n|Intermed|Global'; Val=[1 NaN nc];
        end
        n = spm_input('Null hyp. to assess?','+1','b',But,Val,1);
        if isnan(n)
            if nc == 3,
                n = nc - 1;
            else
                n = nc - spm_input('Effects under null ','0','n1','1',nc-1);
            end
        end
    end
else
    n = 1;
end

%-Enforce orthogonality of multiple contrasts for conjunction
% (Orthogonality within subspace spanned by contrasts)
%--------------------------------------------------------------------------
if nc > 1 && n > 1 && ~spm_FcUtil('|_?',xCon(Ic), xX.xKXs)

    OrthWarn = 0;

    %-Successively orthogonalise
    %-NB: This loop is peculiarly controlled to account for the
    %     possibility that Ic may shrink if some contrasts disappear
    %     on orthogonalisation (i.e. if there are colinearities)
    %----------------------------------------------------------------------
    i = 1;
    while(i < nc), i = i + 1;

        %-Orthogonalise (subspace spanned by) contrast i w.r.t. previous
        %------------------------------------------------------------------
        oxCon = spm_FcUtil('|_',xCon(Ic(i)), xX.xKXs, xCon(Ic(1:i-1)));

        %-See if this orthogonalised contrast has already been entered
        % or is colinear with a previous one. Define a new contrast if
        % neither is the case.
        %------------------------------------------------------------------
        d     = spm_FcUtil('In',oxCon,xX.xKXs,xCon);

        if spm_FcUtil('0|[]',oxCon,xX.xKXs)

            %-Contrast was colinear with a previous one - drop it
            %--------------------------------------------------------------
            Ic(i) = [];
            i     = i - 1;

        elseif any(d)

            %-Contrast unchanged or already defined - note index
            %--------------------------------------------------------------
            Ic(i) = min(d);

        else

            %-Define orthogonalised contrast as new contrast
            %--------------------------------------------------------------
            OrthWarn = OrthWarn + 1;
            conlst     = sprintf('%d,',Ic(1:i-1));
            oxCon.name = sprintf('%s (orth. w.r.t {%s})', xCon(Ic(i)).name,...
                conlst(1:end-1));
            xCon  = [xCon, oxCon];
            Ic(i) = length(xCon);
        end

    end % while...

    if OrthWarn
        warning(sprintf(['Contrasts changed!  %d contrasts orthogonalized ',...
            'to allow conjunction inf.'], OrthWarn))
    end

    SPM.xCon = xCon;
end % if nc>1...
SPM.xCon = xCon;

%-Get contrasts for masking
%--------------------------------------------------------------------------
try
    Im = xSPM.Im;
    if ~isempty(Im)
        maskCon = 1;
    else
        maskCon = 0;
    end
catch
    maskCon = spm_input('mask with other contrast(s)','+1','y/n',[1,0],2);
end
if maskCon
    try
        Im = xSPM.Im;
    catch
        [Im,xCon] = spm_conman(SPM,'T&F',-Inf,...
            'Select contrasts for masking...',' for masking',1);
    end

    %-Threshold for mask (uncorrected p-value)
    %----------------------------------------------------------------------
    try
        pm = xSPM.pm;
    catch
        pm = spm_input('uncorrected mask p-value','+1','r',0.05,1,[0,1]);
    end

    %-Inclusive or exclusive masking
    %----------------------------------------------------------------------
    try
        Ex = xSPM.Ex;
    catch
        Ex = spm_input('nature of mask','+1','b','inclusive|exclusive',[0,1]);
    end
else
    Im = [];
    pm = [];
    Ex = [];
end

%-Create/Get title string for comparison
%--------------------------------------------------------------------------
if nc == 1
    str  = xCon(Ic).name;
else
    str  = [sprintf('contrasts {%d',Ic(1)),sprintf(',%d',Ic(2:end)),'}'];
    if n == nc
        str = [str ' (global null)'];
    elseif n == 1
        str = [str ' (conj. null)'];
    else
        str = [str sprintf(' (Ha: k>=%d)',(nc-n)+1)];
    end
end
if Ex
    mstr = 'masked [excl.] by';
else
    mstr = 'masked [incl.] by';
end
if length(Im) == 1
    str  = sprintf('%s (%s %s at p=%g)',str,mstr,xCon(Im).name,pm);

elseif ~isempty(Im)
    str  = [sprintf('%s (%s {%d',str,mstr,Im(1)),...
        sprintf(',%d',Im(2:end)),...
        sprintf('} at p=%g)',pm)];
end
try
    titlestr = xSPM.title;
    if isempty(titlestr)
        titlestr = str;
    end
catch
    titlestr     = spm_input('title for comparison','+1','s',str);
end

%-Bayesian or classical Inference?
%==========================================================================
if isfield(SPM,'PPM')

    % Make sure SPM.PPM.xCon field exists
    %----------------------------------------------------------------------
    if ~isfield(SPM.PPM,'xCon')
        SPM.PPM.xCon = [];
    end

    % Set Bayesian con type
    %----------------------------------------------------------------------
    if length(SPM.PPM.xCon) < Ic
        SPM.PPM.xCon(Ic).PSTAT = xCon(Ic).STAT;
    end

    % Make this one a Bayesian contrast - if first level
    %----------------------------------------------------------------------
    if isfield(SPM.PPM,'VB')
        [xCon(Ic).STAT] = deal('P');
    end

    if all(strcmp([SPM.PPM.xCon(Ic).PSTAT],'T'))

        % Simple contrast
        %------------------------------------------------------------------
        str = 'Effect size threshold for PPM';

        if isfield(SPM.PPM,'VB') % 1st level Bayes

            % For VB - set default effect size to zero
            %--------------------------------------------------------------
            Gamma = 0;
            xCon(Ic).eidf = spm_input(str,'+1','e',sprintf('%0.2f',Gamma));

        elseif nc == 1 && isempty(xCon(Ic).Vcon) % 2nd level Bayes
            % con image not yet written
            %--------------------------------------------------------------
            if spm_input('Inference',1,'b',{'Bayesian','classical'},[1 0]);

                %-Get Bayesian threshold (Gamma) stored in xCon(Ic).eidf
                % The default is one conditional s.d. of the contrast
                %----------------------------------------------------------
                Gamma         = sqrt(xCon(Ic).c'*SPM.PPM.Cb*xCon(Ic).c);
                xCon(Ic).eidf = spm_input(str,'+1','e',sprintf('%0.2f',Gamma));
                xCon(Ic).STAT = 'P';
            end
        end
    else
        % Compound contrast using Chi^2 statistic
        %------------------------------------------------------------------
        if ~isfield(xCon(Ic),'eidf') || isempty(xCon(Ic).eidf)
            xCon(Ic).eidf = 0; % temporarily
        end
    end
end

%-Compute & store contrast parameters, contrast/ESS images, & SPM images
%==========================================================================
SPM.xCon = xCon;
SPM      = spm_contrasts(SPM, unique([Ic, Im, IcAdd]));
xCon     = SPM.xCon;
STAT     = xCon(Ic(1)).STAT;
VspmSv   = cat(1,xCon(Ic).Vspm);

%-Check conjunctions - Must be same STAT w/ same df
%--------------------------------------------------------------------------
if (nc > 1) && (any(diff(double(cat(1,xCon(Ic).STAT)))) || ...
        any(abs(diff(cat(1,xCon(Ic).eidf))) > 1))
    error('illegal conjunction: can only conjoin SPMs of same STAT & df');
end

%-Degrees of Freedom and STAT string describing marginal distribution
%--------------------------------------------------------------------------
df     = [xCon(Ic(1)).eidf xX.erdf];
if nc > 1
    if n > 1
        str = sprintf('^{%d \\{Ha:k\\geq%d\\}}',nc,(nc-n)+1);
    else
        str = sprintf('^{%d \\{Ha:k=%d\\}}',nc,(nc-n)+1);
    end
else
    str = '';
end

switch STAT
    case 'T'
        STATstr = sprintf('%c%s_{%.0f}','T',str,df(2));
    case 'F'
        STATstr = sprintf('%c%s_{%.0f,%.0f}','F',str,df(1),df(2));
    case 'P'
        STATstr = sprintf('%s^{%0.2f}','PPM',df(1));
end

%-Compute (unfiltered) SPM pointlist for masked conjunction requested
%==========================================================================
fprintf('\t%-32s: %30s','SPM computation','...initialising')            %-#

%-Compute conjunction as minimum of SPMs
%--------------------------------------------------------------------------
Z     = Inf;
for i = Ic
    Z = min(Z,spm_get_data(xCon(i).Vspm,XYZ));
end

%-Copy of Z and XYZ before masking, for later use with FDR
%--------------------------------------------------------------------------
XYZum = XYZ;
Zum   = Z;

%-Compute mask and eliminate masked voxels
%--------------------------------------------------------------------------
for i = Im
    fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...masking')           %-#

    Mask = spm_get_data(xCon(i).Vspm,XYZ);
    um   = spm_u(pm,[xCon(i).eidf,xX.erdf],xCon(i).STAT);
    if Ex
        Q = Mask <= um;
    else
        Q = Mask >  um;
    end
    XYZ   = XYZ(:,Q);
    Z     = Z(Q);
    if isempty(Q)
        fprintf('\n')                                                   %-#
        warning(sprintf('No voxels survive masking at p=%4.2f',pm));
        break
    end
end

%==========================================================================
% - H E I G H T   &   E X T E N T   T H R E S H O L D S
%==========================================================================

u   = -Inf;        % height threshold
k   = 0;           % extent threshold {voxels}

%-Get FDR mode
%--------------------------------------------------------------------------
try
    topoFDR = spm_get_defaults('stats.topoFDR');
catch
    topoFDR = true;
end

%-Height threshold - classical inference
%--------------------------------------------------------------------------
if STAT ~= 'P'

    %-Get height threshold
    %----------------------------------------------------------------------
    fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...height threshold')  %-#
    try
        thresDesc = xSPM.thresDesc;
    catch
        if topoFDR
            str = 'FWE|none';
        else
            str = 'FWE|FDR|none';
        end
        thresDesc = spm_input('p value adjustment to control','+1','b',str,[],1);
    end

    switch thresDesc

        case 'FWE' % Family-wise false positive rate
            %--------------------------------------------------------------
            try
                u = xSPM.u;
            catch
                u = spm_input('p value (FWE)','+0','r',0.05,1,[0,1]);
            end
            thresDesc = ['p<' num2str(u) ' (' thresDesc ')'];
            u = spm_uc(u,df,STAT,R,n,S);

        case 'FDR' % False discovery rate
            %--------------------------------------------------------------
            if topoFDR,
                fprintf('\n');                                          %-#
                error('Change defaults.stats.topoFDR to use voxel FDR.');
            end
            try
                u = xSPM.u;
            catch
                u = spm_input('p value (FDR)','+0','r',0.05,1,[0,1]);
            end
            thresDesc = ['p<' num2str(u) ' (' thresDesc ')'];
            u = spm_uc_FDR(u,df,STAT,n,VspmSv,0);

        case 'none'  % No adjustment: p for conjunctions is p of the conjunction SPM
            %--------------------------------------------------------------
            try
                u = xSPM.u;
            catch
                u = spm_input(['threshold {',STAT,' or p value}'],'+0','r',0.001,1);
            end
            if u <= 1
                thresDesc = ['p<' num2str(u) ' (unc.)'];
                u = spm_u(u^(1/n),df,STAT);
            else
                thresDesc = [STAT '=' num2str(u) ];
            end

        otherwise
            %--------------------------------------------------------------
            fprintf('\n');                                              %-#
            error(sprintf('Unknown control method "%s".',thresDesc));

    end % switch thresDesc

    %-Compute p-values for topological and voxel-wise FDR (all search voxels)
    %----------------------------------------------------------------------
    if ~topoFDR
        %-Voxel-wise FDR
        fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...for voxelFDR')  %-#
        switch STAT
            case 'Z'
                Ps   = (1-spm_Ncdf(Zum)).^n;
            case 'T'
                Ps   = (1 - spm_Tcdf(Zum,df(2))).^n;
            case 'X'
                Ps   = (1-spm_Xcdf(Zum,df(2))).^n;
            case 'F'
                Ps   = (1 - spm_Fcdf(Zum,df)).^n;
        end
        Ps = sort(Ps);
    end

    %-Peak FDR
    %----------------------------------------------------------------------
    [up,Pp] = spm_uc_peakFDR(0.05,df,STAT,R,n,Zum,XYZum,u);

    %-Cluster FDR
    %----------------------------------------------------------------------
    if STAT == 'T' && n == 1
        V2R        = 1/prod(SPM.xVol.FWHM(SPM.xVol.DIM > 1));
        [uc,Pc,ue] = spm_uc_clusterFDR(0.05,df,STAT,R,n,Zum,XYZum,V2R,u);
    else
        uc  = NaN;
        ue  = NaN;
        Pc  = [];
    end

    %-Peak FWE
    %----------------------------------------------------------------------
    uu      = spm_uc(0.05,df,STAT,R,n,S);

%-Height threshold - Bayesian inference
%--------------------------------------------------------------------------
elseif STAT == 'P'

    u_default = 1 - 1/SPM.xVol.S;
    str       = 'Posterior probability threshold for PPM';
    u         = spm_input(str,'+0','r',u_default,1);
    thresDesc = ['P>'  num2str(u) ' (PPM)'];

end % (if STAT)

%-Calculate height threshold filtering
%--------------------------------------------------------------------------
Q      = find(Z > u);

%-Apply height threshold
%--------------------------------------------------------------------------
Z      = Z(:,Q);
XYZ    = XYZ(:,Q);
if isempty(Q)
    fprintf('\n');                                                      %-#
    warning(sprintf('No voxels survive height threshold u=%0.2g',u))
end

%-Extent threshold (disallowed for conjunctions)
%--------------------------------------------------------------------------
if ~isempty(XYZ) && nc == 1

    fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...extent threshold')  %-#

    %-Get extent threshold [default = 0]
    %----------------------------------------------------------------------
    try
        k = xSPM.k;
    catch
        k = spm_input('& extent threshold {voxels}','+1','r',0,1,[0,Inf]);
    end

    %-Calculate extent threshold filtering
    %----------------------------------------------------------------------
    A     = spm_clusters(XYZ);
    Q     = [];
    for i = 1:max(A)
        j = find(A == i);
        if length(j) >= k; Q = [Q j]; end
    end

    % ...eliminate voxels
    %----------------------------------------------------------------------
    Z     = Z(:,Q);
    XYZ   = XYZ(:,Q);
    if isempty(Q)
        fprintf('\n');                                                  %-#
        warning(sprintf('No voxels survive extent threshold k=%0.2g',k))
    end

else

    k = 0;

end % (if ~isempty(XYZ))

%-For Bayesian inference provide (default) option to display contrast values
%--------------------------------------------------------------------------
if STAT == 'P'
    if spm_input('Plot effect-size/statistic',1,'b',{'Yes','No'},[1 0])
        Z = spm_get_data(xCon(Ic).Vcon,XYZ);
    end
end

%==========================================================================
% - E N D
%==========================================================================
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done')                %-#
spm('Pointer','Arrow')

%-Assemble output structures of unfiltered data
%==========================================================================
xSPM   = struct( ...
            'swd',      swd,...
            'title',    titlestr,...
            'Z',        Z,...
            'n',        n,...
            'STAT',     STAT,...
            'df',       df,...
            'STATstr',  STATstr,...
            'Ic',       Ic,...
            'Im',       Im,...
            'pm',       pm,...
            'Ex',       Ex,...
            'u',        u,...
            'k',        k,...
            'XYZ',      XYZ,...
            'XYZmm',    SPM.xVol.M(1:3,:)*[XYZ; ones(1,size(XYZ,2))],...
            'S',        SPM.xVol.S,...
            'R',        SPM.xVol.R,...
            'FWHM',     SPM.xVol.FWHM,...
            'M',        SPM.xVol.M,...
            'iM',       SPM.xVol.iM,...
            'DIM',      SPM.xVol.DIM,...
            'VOX',      VOX,...
            'Vspm',     VspmSv,...
            'thresDesc',thresDesc);

% RESELS per voxel (density) if it exists
%--------------------------------------------------------------------------
xSPM.VRpv = SPM.xVol.VRpv;
try
    xSPM.units = SPM.xVol.units;
catch
    try, xSPM.units = varargin{1}.units; end;
end

% p-values for topological and voxel-wise FDR
%--------------------------------------------------------------------------
try, xSPM.Ps    = Ps;             end  % voxel   FDR
try, xSPM.Pp    = Pp;             end  % peak    FDR
try, xSPM.Pc    = Pc;             end  % cluster FDR

% 0.05 critical thresholds for FWEp, FDRp, FWEc, FDRc
%--------------------------------------------------------------------------
try, xSPM.uc    = [uu up ue uc];  end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function fx = circle ( n, x )

%*****************************************************************************80
%
%% CIRCLE evaluates the circle function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N,2), the point coordinates.
%
%    Output, real FX(N), the function values.
%
  global c
  global r

  if ( isempty ( c ) )
    c = [ 0.0, 0.5 ];
  end

  if ( isempty ( r ) )
    r = 0.5;
  end

  nc = ones ( n, 1 ) * c;
  nr = ones ( n, 1 ) * r;

  fx = nr - sqrt ( sum ( ( x - nc ).^2, 2 ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function title = p40_title ( )

%*****************************************************************************80
%
%% P40_TITLE returns a title for problem 40.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The Bohachevsky Function #3';

  return
end


text:
classdef BRISK < handle
    %BRISK  Class implementing the BRISK keypoint detector and descriptor extractor
    %
    % As described in [LCS11].
    %
    % ## References
    % [LCS11]:
    % > Stefan Leutenegger, Margarita Chli, and Roland Yves Siegwart.
    % > "BRISK: Binary Robust Invariant Scalable Keypoints".
    % > In IEEE International Conference on Computer Vision (ICCV) 2011,
    % > pages 2548-2555. IEEE, 2011.
    %
    % See also: cv.BRISK.BRISK, cv.FeatureDetector, cv.DescriptorExtractor,
    %  detectBRISKFeatures, extractFeatures, BRISKPoints
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = BRISK(varargin)
            %BRISK  The BRISK constructor
            %
            %     obj = cv.BRISK()
            %     obj = cv.BRISK(radiusList, numberList)
            %     obj = cv.BRISK(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __radiusList__ defines the radii (in pixels) where the samples
            %   around a keypoint are taken (for keypoint scale 1).
            % * __numberList__ defines the number of sampling points on the
            %   sampling circle. Must be the same size as `radiusList`.
            %
            % ## Options
            % Options accepted by first variant:
            %
            % * __Threshold__ FAST/AGAST detection threshold score. default 30
            % * __Octaves__ detection octaves. Use 0 to do single scale.
            %   default 3
            % * __PatternScale__ apply this scale to the pattern used for
            %   sampling the neighbourhood of a keypoint. default 1.0
            %
            % Options accepted by second variant for a custom pattern:
            %
            % * __Threshold__ FAST/AGAST detection threshold score. default 20
            % * __Octaves__ detection octaves. Use 0 to do single scale.
            %   default 3
            % * __DMax__ threshold for the short pairings used for descriptor
            %   formation (in pixels for keypoint scale 1). default 5.85
            % * __DMin__ threshold for the long pairings used for orientation
            %   determination (in pixels for keypoint scale 1). default 8.2
            % * __IndexChange__ index remapping of the bits.
            %   default empty vector.
            %
            % See also: cv.BRISK.detectAndCompute
            %
            this.id = BRISK_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.BRISK
            %
            if isempty(this.id), return; end
            BRISK_(this.id, 'delete');
        end

        function typename = typeid(this)
            %TYPEID  Name of the C++ type (RTTI)
            %
            %     typename = obj.typeid()
            %
            % ## Output
            % * __typename__ Name of C++ type
            %
            typename = BRISK_(this.id, 'typeid');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.BRISK.empty, cv.BRISK.load
            %
            BRISK_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.BRISK.clear, cv.BRISK.load
            %
            b = BRISK_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.BRISK.load
            %
            BRISK_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.BRISK.save
            %
            BRISK_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.BRISK.save, cv.BRISK.load
            %
            name = BRISK_(this.id, 'getDefaultName');
        end
    end

    %% Features2D: FeatureDetector + DescriptorExtractor
    methods
        function ntype = defaultNorm(this)
            %DEFAULTNORM  Returns the default norm type
            %
            %     ntype = obj.defaultNorm()
            %
            % ## Output
            % * __ntype__ Norm type. One of `cv::NormTypes`:
            %   * __Inf__
            %   * __L1__
            %   * __L2__
            %   * __L2Sqr__
            %   * __Hamming__
            %   * __Hamming2__
            %
            % Always `Hamming` for BRISK.
            %
            % See also: cv.BRISK.compute, cv.DescriptorMatcher
            %
            ntype = BRISK_(this.id, 'defaultNorm');
        end

        function sz = descriptorSize(this)
            %DESCRIPTORSIZE  Returns the descriptor size in bytes
            %
            %     sz = obj.descriptorSize()
            %
            % ## Output
            % * __sz__ Descriptor size.
            %
            % See also: cv.BRISK.descriptorType, cv.BRISK.compute
            %
            sz = BRISK_(this.id, 'descriptorSize');
        end

        function dtype = descriptorType(this)
            %DESCRIPTORTYPE  Returns the descriptor type
            %
            %     dtype = obj.descriptorType()
            %
            % ## Output
            % * __dtype__ Descriptor type, one of numeric MATLAB class names.
            %
            % Always `uint8` for BRISK.
            %
            % See also: cv.BRISK.descriptorSize, cv.BRISK.compute
            %
            dtype = BRISK_(this.id, 'descriptorType');
        end

        function keypoints = detect(this, img, varargin)
            %DETECT  Detects keypoints in an image or image set
            %
            %     keypoints = obj.detect(img)
            %     keypoints = obj.detect(imgs)
            %     [...] = obj.detect(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ Image (first variant), 8-bit grayscale image.
            % * __imgs__ Image set (second variant), cell array of images.
            %
            % ## Output
            % * __keypoints__ The detected keypoints. In the first variant, a
            %   1-by-N structure array. In the second variant of the method,
            %   `keypoints{i}` is a set of keypoints detected in `imgs{i}`.
            %
            % ## Options
            % * __Mask__ A mask specifying where to look for keypoints
            %   (optional). It must be a logical or 8-bit integer matrix with
            %   non-zero values in the region of interest. In the second
            %   variant, it is a cell-array of masks for each input image,
            %   `masks{i}` is a mask for `imgs{i}`. Not set by default.
            %
            % See also: cv.BRISK.compute, cv.BRISK.detectAndCompute
            %
            keypoints = BRISK_(this.id, 'detect', img, varargin{:});
        end

        function [descriptors, keypoints] = compute(this, img, keypoints)
            %COMPUTE  Computes the descriptors for a set of keypoints detected in an image or image set
            %
            %     [descriptors, keypoints] = obj.compute(img, keypoints)
            %     [descriptors, keypoints] = obj.compute(imgs, keypoints)
            %
            % ## Input
            % * __img__ Image (first variant), 8-bit grayscale image.
            % * __imgs__ Image set (second variant), cell array of images.
            % * __keypoints__ Input collection of keypoints. Keypoints for
            %   which a descriptor cannot be computed are removed. Sometimes
            %   new keypoints can be added, for example: cv.SIFT duplicates
            %   keypoint with several dominant orientations (for each
            %   orientation). In the first variant, this is a struct-array of
            %   detected keypoints. In the second variant, it is a cell-array,
            %   where `keypoints{i}` is a set of keypoints detected in
            %   `imgs{i}` (a struct-array like before).
            %
            % ## Output
            % * __descriptors__ Computed descriptors. In the second variant of
            %   the method `descriptors{i}` are descriptors computed for a
            %   `keypoints{i}`. Row `j` in `descriptors` (or `descriptors{i}`)
            %   is the descriptor for `j`-th keypoint.
            % * __keypoints__ Optional output with possibly updated keypoints.
            %
            % See also: cv.BRISK.detect, cv.BRISK.detectAndCompute
            %
            [descriptors, keypoints] = BRISK_(this.id, 'compute', img, keypoints);
        end

        function [keypoints, descriptors] = detectAndCompute(this, img, varargin)
            %DETECTANDCOMPUTE  Detects keypoints and computes their descriptors
            %
            %     [keypoints, descriptors] = obj.detectAndCompute(img)
            %     [...] = obj.detectAndCompute(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ Image, input 8-bit grayscale image.
            %
            % ## Output
            % * __keypoints__ The detected keypoints. A 1-by-N structure array
            %   with the following fields:
            %   * __pt__ coordinates of the keypoint `[x,y]`
            %   * __size__ diameter of the meaningful keypoint neighborhood
            %   * __angle__ computed orientation of the keypoint (-1 if not
            %     applicable); it's in [0,360) degrees and measured relative
            %     to image coordinate system (y-axis is directed downward),
            %     i.e in clockwise.
            %   * __response__ the response by which the most strong keypoints
            %     have been selected. Can be used for further sorting or
            %     subsampling.
            %   * __octave__ octave (pyramid layer) from which the keypoint
            %     has been extracted.
            %   * **class_id** object class (if the keypoints need to be
            %     clustered by an object they belong to).
            % * __descriptors__ Computed descriptors. Output concatenated
            %   vectors of descriptors. Each descriptor is a vector of length
            %   cv.BRISK.descriptorSize, so the total size of descriptors will
            %   be `numel(keypoints) * obj.descriptorSize()`. A matrix of size
            %   N-by-sz of class `uint8`, one row per keypoint.
            %
            % ## Options
            % * __Mask__ optional mask specifying where to look for keypoints.
            %   Not set by default.
            % * __Keypoints__ If passed, then the method will use the provided
            %   vector of keypoints instead of detecting them, and the
            %   algorithm just computes their descriptors.
            %
            % See also: cv.BRISK.detect, cv.BRISK.compute
            %
            [keypoints, descriptors] = BRISK_(this.id, 'detectAndCompute', img, varargin{:});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function pdf = beta_binomial_pdf ( x, a, b, c )

%*****************************************************************************80
%
%% BETA_BINOMIAL_PDF evaluates the Beta Binomial PDF.
%
%  Discussion:
%
%    PDF(X)(A,B,C) = Beta(A+X,B+C-X)
%      / ( (C+1) * Beta(X+1,C-X+1) * Beta(A,B) )  for 0 <= X <= C.
%
%    This PDF can be reformulated as:
%
%      The beta binomial probability density function for X successes
%      out of N trials is
%
%      PDF2(X)( N, MU, THETA ) =
%        C(N,X) * Product ( 0 <= R <= X - 1 ) ( MU + R * THETA )
%               * Product ( 0 <= R <= N - X - 1 ) ( 1 - MU + R * THETA )
%               / Product ( 0 <= R <= N - 1 )  ( 1 + R * THETA )
%
%      where
%
%        C(N,X) is the combinatorial coefficient;
%        MU is the expectation of the underlying Beta distribution;
%        THETA is a shape parameter.
%
%      A THETA value of 0 ( or A+B --> Infinity ) results in the binomial
%      distribution:
%
%        PDF2(X) ( N, MU, 0 ) = C(N,X) * MU**X * ( 1 - MU )**(N-X)
%
%    Given A, B, C for PDF, then the equivalent PDF2 has:
%
%      N     = C
%      MU    = A / ( A + B )
%      THETA = 1 / ( A + B )
%
%    Given N, MU, THETA for PDF2, the equivalent PDF has:
%
%      A = MU / THETA
%      B = ( 1 - MU ) / THETA
%      C = N
%
%    BETA_BINOMIAL_PDF(X)(1,1,C) = UNIFORM_DISCRETE_PDF(X)(0,C-1)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, the argument of the PDF.
%
%    Input, real A, B, parameters of the PDF.
%    0.0D+00 < A,
%    0.0D+00 < B.
%
%    Input, integer C, a parameter of the PDF.
%    0 <= C.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < 0 )

    pdf = 0.0;

  elseif ( x <= c )

    pdf = beta ( a + x, b + c - x ) ...
      / ( ( c + 1 ) * beta ( x + 1, c - x + 1 ) * beta ( a, b ) );

  elseif ( c < x )

    pdf = 0.0;

  end

  return
end


text:
function [ x, indx, k, nstack, stacks, ncan ] = r8vec_backtrack ( n, ...
  x, indx, k, nstack, stacks, ncan )

%*****************************************************************************80
%
%% R8VEC_BACKTRACK supervises a backtrack search for an R8VEC.
%
%  Discussion:
%
%    The routine tries to construct a real vector one index at a time,
%    using possible candidates as supplied by the user.
%
%    At any time, the partially constructed vector may be discovered to be
%    unsatisfactory, but the routine records information about where the
%    last arbitrary choice was made, so that the search can be
%    carried out efficiently, rather than starting out all over again.
%
%    First, call the routine with INDX = 0 so it can initialize itself.
%
%    Now, on each return from the routine, if INDX is:
%      1, you've just been handed a complete candidate vector;
%         Admire it, analyze it, do what you like.
%      2, please determine suitable candidates for position X(K).
%         Return the number of candidates in NCAN(K), adding each
%         candidate to the end of STACKS, and increasing NSTACK.
%      3, you're done.  Stop calling the routine;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2011
%
%  Author:
%
%    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf.
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the number of positions to be filled in the vector.
%
%    Input, real X(N), the partially filled in candidate vector.
%
%    Input, integer INDX, a communication flag.
%    * 0, to begin a backtracking search.
%    * 2, the requested candidates for position K have been added to
%      STACKS, and NCAN(K) was updated.
%
%    Input, integer K, the index in X that we are trying to fill.
%
%    Input, integer NSTACK, the current length of the stack.
%
%    Input, real STACKS(MAXSTACK), a list of all current candidates for
%    all positions 1 through K.
%
%    Input, integer NCAN(N), lists the current number of candidates for
%    all positions 1 through K.
%
%    Output, real X(N), the partially filled in candidate vector.
%
%    Output, integer INDX, a communication flag.
%    * 1, a complete output vector has been determined and returned in X(1:N);
%    * 2, candidates are needed for position X(K);
%    * 3, no more possible vectors exist.
%
%    Output, integer K, the index in X that we are trying to fill.
%
%    Output, integer NSTACK, the current length of the stack.
%
%    Output, real STACKS(MAXSTACK), a list of all current candidates for
%    all positions 1 through K.
%
%    Output, integer NCAN(N), lists the current number of candidates for
%    all positions 1 through K.
%

%
%  If this is the first call, request a candidate for position 1.
%
  if ( indx == 0 )
    k = 1;
    nstack = 0;
    indx = 2;
    return
  end
%
%  Examine the stack.
%
  while ( 1 )
%
%  If there are candidates for position K, take the first available
%  one off the stack, and increment K.
%
%  This may cause K to reach the desired value of N, in which case
%  we need to signal the user that a complete set of candidates
%  is being returned.
%
    if ( 0 < ncan(k) )

      x(k) = stacks(nstack);
      nstack = nstack - 1;
      ncan(k) = ncan(k) - 1;

      if ( k ~= n )
        k = k + 1;
        indx = 2;
      else
        indx = 1;
      end

      break
%
%  If there are no candidates for position K, then decrement K.
%  If K is still positive, repeat the examination of the stack.
%
    else

      k = k - 1;

      if ( k <= 0 )
        indx = 3;
        break
      end

    end

  end

  return
end


text:
% function [Lh,Ph,Mu,Pi,LL]=mfa(X,M,K,cyc,tol);
%
% Maximum Likelihood Mixture of Factor Analysis using EM
%
% X - data matrix
% M - number of mixtures (default 1)
% K - number of factors in each mixture (default 2)
% cyc - maximum number of cycles of EM (default 100)
% tol - termination tolerance (prop change in likelihood) (default 0.0001)
%
% Lh - factor loadings
% Ph - diagonal uniquenesses matrix
% Mu - mean vectors
% Pi - priors
% LL - log likelihood curve
%
% Iterates until a proportional change < tol in the log likelihood
% or cyc steps of EM

function [Lh, Ph,  Mu, Pi, LL] = mfa(X,M,K,cyc,tol)

if nargin<5   tol=0.0001; end;
if nargin<4   cyc=100; end;
if nargin<3   K=2; end;
if nargin<2   M=1; end;

N=length(X(:,1));
D=length(X(1,:));
tiny=exp(-700);

%rand('state',0);

fprintf('\n');

if (M==1)
  [Lh,Ph,LL]=ffa(X,K,cyc,tol);
  Mu=mean(X);
  Pi=1;
else
  if N==1
    mX = X;
  else
    mX=mean(X);
  end
  cX=cov(X);
  scale=det(cX)^(1/D);
  randn('state',0);
  Lh=randn(D*M,K)*sqrt(scale/K);
  Ph=diag(cX)+tiny;
  Pi=ones(M,1)/M;
  %randn('state',0);
  Mu=randn(M,D)*sqrtm(cX)+ones(M,1)*mX;
  oldMu=Mu;
  I=eye(K);

  lik=0;
  LL=[];

  H=zeros(N,M); 	% E(w|x)
  EZ=zeros(N*M,K);
  EZZ=zeros(K*M,K);
  XX=zeros(D*M,D);
  s=zeros(M,1);
  const=(2*pi)^(-D/2);
  %%%%%%%%%%%%%%%%%%%%
  for i=1:cyc;

    %%%% E Step %%%%

    Phi=1./Ph;
    Phid=diag(Phi);
    for k=1:M
      Lht=Lh((k-1)*D+1:k*D,:);
      LP=Phid*Lht;
      MM=Phid-LP*inv(I+Lht'*LP)*LP';
      dM=sqrt(det(MM));
      Xk=(X-ones(N,1)*Mu(k,:));
      XM=Xk*MM;
      H(:,k)=const*Pi(k)*dM*exp(-0.5*rsum(XM.*Xk));
      EZ((k-1)*N+1:k*N,:)=XM*Lht;
    end;

    Hsum=rsum(H);
    oldlik=lik;
    lik=sum(log(Hsum+(Hsum==0)*exp(-744)));

    Hzero=(Hsum==0); Nz=sum(Hzero);
    H(Hzero,:)=tiny*ones(Nz,M)/M;
    Hsum(Hzero)=tiny*ones(Nz,1);

    H=rdiv(H,Hsum);
    s=csum(H);
    s=s+(s==0)*tiny;
    s2=sum(s)+tiny;

    for k=1:M
      kD=(k-1)*D+1:k*D;
      Lht=Lh(kD,:);
      LP=Phid*Lht;
      MM=Phid-LP*inv(I+Lht'*LP)*LP';
      Xk=(X-ones(N,1)*Mu(k,:));
      XX(kD,:)=rprod(Xk,H(:,k))'*Xk/s(k);
      beta=Lht'*MM;
      EZZ((k-1)*K+1:k*K,:)=I-beta*Lht +beta*XX(kD,:)*beta';
    end;

    %%%% log likelihood %%%%

    LL=[LL lik];
    fprintf('cycle %g   \tlog likelihood %g ',i,lik);

    if (i<=2)
      likbase=lik;
    elseif (lik<oldlik)
      fprintf(' violation');
    elseif ((lik-likbase)<(1 + tol)*(oldlik-likbase)||~isfinite(lik))
      break;
    end;

    fprintf('\n');

    %%%% M Step %%%%

    % means and covariance structure

    Ph=zeros(D,1);
    for k=1:M
      kD=(k-1)*D+1:k*D;
      kK=(k-1)*K+1:k*K;
      kN=(k-1)*N+1:k*N;

      T0=rprod(X,H(:,k));
      T1=T0'*[EZ(kN,:) ones(N,1)];
      XH=EZ(kN,:)'*H(:,k);
      T2=inv([s(k)*EZZ(kK,:) XH; XH' s(k)]);
      T3=T1*T2;
      Lh(kD,:)=T3(:,1:K);
      Mu(k,:)=T3(:,K+1)';
      T4=diag(T0'*X-T3*T1')/s2;
      Ph=Ph+T4.*(T4>0);
    end;

    Phmin=exp(-700);
    Ph=Ph.*(Ph>Phmin)+(Ph<=Phmin)*Phmin; % to avoid zero variances

    % priors
    Pi=s'/s2;

  end;
  fprintf('\n');
end;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function Y = fillnans(varargin)
% FILLNANS replaces all NaNs in array using inverse-distance weighting.
%
% Y = FILLNANS(X) replaces all NaNs in the vector or array X by
% inverse-distance weighted interpolation:
%                       Y = sum(X/D^3)/sum(1/D^3)
% where D is the distance (in pixels) from the NaN node to all non-NaN
% values X. Values farther from a known non-NaN value will tend toward the
% average of all the values.
%
% Y = FILLNANS(...,'power',p) uses a power of p in the weighting
% function. The higer the value of p, the stronger the weighting.
%
% Y = FILLNANS(...,'radius',d) only used pixels < d pixels away in
% for weighted averaging.
%
% NOTE: Use in conjunction with INVDISTGRID to grid and interpolate x,y,z
% data.
%
% See also INPAINT_NANS
%
% Ian M. Howat, Applied Physics Lab, University of Washington
% ihowat@apl.washington.edu Ian M. Howat
% Version 1: 05-Jul-2007 17:28:57
%   Revision 1: 16-Jul-2007 17:47:43
%       Added increment expression to waitbar to reduce number of times its
%       called. Provided by John D'Errico.
%   Revision 2: 16-Jul-2007 18:40:34
%       Added radius option and 'option',value varargin parser.
%   Revision 3: 18-Jul-2007 10:25:23
%       Adopted several code efficiency revisions made by Urs, including
%       removing the waitbar.

%parse input and set defualts:
X = varargin{1}; %input array
Y = X; %output array
n = 2; %weighting power
d = 0; %distance cut-off radius (0= all pixels, no cut-off)
if nargin > 1 && nargin < 6
    for k=2:2:length(varargin)
        if isnumeric(varargin{k}) || ~isnumeric(varargin{k+1})
            error('Input arguments must be in ''option'',value form.')
        end
        switch lower(varargin{k})	% (Urs:less error prone)
            case 'power'
                n = varargin{k+1};
            case 'radius'
                d = varargin{k+1};
            otherwise
                error(['Unrecognized input argument: ',varargin{k}])
        end
    end
elseif nargin >= 6
    error('Too many input arguments')
end

%(Urs:use ISNAN()/NOT() once only)
ix=isnan(X);
[rn,cn]=find(ix);   %row,col of nans
ix=~ix;
[r,c]=find(ix);     %row,col of non-nans
ind=find(ix);       %index of non-nans

%Break distance-finding loops into with cut-off and without cut-off
%versions. The cutoff conditional statement adds time
%if cut-off values near the max pixel distance are used.

if d %distance cut-off loop
    d=d.^2;				% (Urs:allows first step without SQRT())
    for k = 1:length(rn)
        D = (rn(k)-r).^2+(cn(k)-c).^2;	% (Urs:no SQRT() here)
        Dd = D < d;
        if sum(Dd) ~= 0
            D=1./sqrt(D(Dd)).^n; %(Urs: Compute once only for valid <D>s)
            Y(rn(k),cn(k)) = sum(X(ind(Dd)).*D)./ sum(D);
        end
    end
else %no distance cut-off loop
    for k = 1:length(rn)
        D = 1./(sqrt((rn(k)-r).^2+(cn(k)-c).^2)).^n;% (Urs:compute once only)
        Y(rn(k),cn(k)) = sum(X(ind).*D)./sum(D);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function y=qrpermute(x,p)
%QRPERMUTE transpose or permute a quaternion array y=[x,p]
%
% Inputs:   x(4m,...)  Real quaternions array
%           p          new order of dimensions [default [2 1 ...]
%
% Outputs:  y(4n,...)  output real quaternion array

%      Copyright (C) Mike Brookes 2012
%      Version: $Id: qrpermute.m 1619 2012-03-15 09:31:31Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s=size(x);
if nargin<2
    p=[2 1 3:length(s)];
end
s(1)=s(1)/4;
t=s(p);
t(1)=4*t(1);
y=reshape(permute(reshape(x,[4 s]),[1 p+1]),t);

author: jtkim-kaist

repo: VAD

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


text:
function [ v, seed ] = vibonacci ( n, seed )

%*****************************************************************************80
%
%% VIBONACCI computes the first N Vibonacci numbers.
%
%  Discussion:
%
%    The "Vibonacci numbers" are a generalization of the Fibonacci numbers:
%      V(N+1) = +/- V(N) +/- V(N-1)
%    where the signs are chosen randomly.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Brian Hayes,
%    The Vibonacci Numbers,
%    American Scientist,
%    July-August 1999, Volume 87, Number 4.
%
%    Divakar Viswanath,
%    Random Fibonacci sequences and the number 1.13198824,
%    Mathematics of Computation, 1998.
%
%  Parameters:
%
%    Input, integer N, the highest number to compute.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, integer V(N), the first N Vibonacci numbers.  By convention,
%    V(1) and V(2) are taken to be 1.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  if ( n <= 0 )
    v = [];
    return
  end

  v(1) = 1;

  if ( n <= 1 )
    return
  end

  v(2) = 1;

  for i = 3 : n

    [ j, seed ] = i4_uniform_ab ( 0, 1, seed );

    if ( j == 0 )
      s1 = -1;
    else
      s1 = +1;
    end

    [ j, seed ] = i4_uniform_ab ( 0, 1, seed );

    if ( j == 0 )
      s2 = -1;
    else
      s2 = +1;
    end

    v(i) = s1 * v(i-1) + s2 * v(i-2);

  end

  return
end


text:
function [centroid, dimension_weight, class] = Entropy_Weighting_Subspace_Kmeans(data, iteration, K, beta, lambda, verbose)

% Pre-allocate weights
centroid = zeros(K,size(data,2)); % centroid for each class
class = zeros(size(data,1),1); % classification result for each observation
dimension_weight = ones(K,size(data,2)); % weights for each dimension in each class

% Initialization
n = size(data,1);
J = 0; % used to record objective function value

% First round
centroid = data(unidrnd(n,K,1),:); % randomly choose K initial centroids
dimension_weight = 1/size(data,2) * dimension_weight; % initialize weights for each dimension using uniform distribution

% Partially optimization
for i = 1:iteration

    J0 = J;

    % body part
    class = classify(centroid,dimension_weight,data);
    centroid = centroid_update(centroid,data,class);
    [dimension_weight,J] = dimension_weight_update(K,centroid,class,dimension_weight,data,beta,lambda);

    % whether display intermediate information
    if(verbose == 1)
        disp(['Objective Function(J) Value: ',num2str(J)]);
    end

    % early stop condition
    if(abs(J-J0)<1e-9)
        fprintf('*** Clustering terminates after %i iterations ***\n',i);
        break;
    end
end
end

% Details on the scheme of updating dimension weights
function [alpha,J] = dimension_weight_update(K,m,c,alpha,data,beta)

  % Under writting...

end

function [result] = classify(m,alpha,data)
result = zeros(size(data,1),1);
% Construct temporary matrix for efficiently computing dimention-weighted distance
matrix = zeros(size(m,1),size(data,2));
for i = 1:size(result,1)
    for j = 1:size(m,1)
        matrix(j,:) = (m(j,:)-data(i,:)) .* (m(j,:)-data(i,:));
    end
    temp = sum(alpha .* matrix,2);
    % To avoid more than one class having the minimum distance.
    t_index = find(temp == min(temp));
    result(i,1) = t_index(1,1);
end
end

function [m] = centroid_update(m,data,c)
for i = 1:size(m,1)
    if(~isempty(data(c==i,:)))
        m(i,:) = mean(data(c==i,:),1);
    else
        continue;
    end
end
end

author: xuyxu

repo: Clustering

save_path: github-repos/MATLAB/xuyxu-Clustering


text:
function [y,r,vr]=ssa(x1,L)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% -----------------------------------------------------------------
%    Author: Francisco Javier Alonso Sanchez    e-mail:fjas@unex.es
%    Departament of Electronics and Electromecanical Engineering
%    Industrial Engineering School
%    University of Extremadura
%    Badajoz
%    Spain
% -----------------------------------------------------------------
%
% SSA generates a trayectory matrix X from the original series x1
% by sliding a window of length L. The trayectory matrix is aproximated
% using Singular Value Decomposition. The last step reconstructs
% the series from the aproximated trayectory matrix. The SSA applications
% include smoothing, filtering, and trend extraction.
% The algorithm used is described in detail in: Golyandina, N., Nekrutkin,
% V., Zhigljavsky, A., 2001. Analisys of Time Series Structure - SSA and
% Related Techniques. Chapman & Hall/CR.

% x1 Original time series (column vector form)
% L  Window length
% y  Reconstructed time series
% r  Residual time series r=x1-y
% vr Relative value of the norm of the approximated trajectory matrix with respect
%	  to the original trajectory matrix

% The program output is the Singular Spectrum of x1 (must be a column vector),
% using a window length L. You must choose the components be used to reconstruct
%the series in the form [i1,i2:ik,...,iL], based on the Singular Spectrum appearance.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Step1 : Build trayectory matrix

   N=length(x1);
   if L>N/2;L=N-L;end
	K=N-L+1;
   X=zeros(L,K);
	for i=1:K
	  X(1:L,i)=x1(i:L+i-1);
	end

% Step 2: SVD

   S=X*X';
	[U,autoval]=eig(S);
	[d,i]=sort(-diag(autoval));
   d=-d;
   U=U(:,i);sev=sum(d);
	plot((d./sev)*100),hold on,plot((d./sev)*100,'rx');
	title('Singular Spectrum');xlabel('Eigenvalue Number');ylabel('Eigenvalue (% Norm of trajectory matrix retained)')
   V=(X')*U;
   rc=U*V';

% Step 3: Grouping

   I=input('Choose the agrupation of components to reconstruct the series in the form I=[i1,i2:ik,...,iL]  ')
   Vt=V';
   rca=U(:,I)*Vt(I,:);

% Step 4: Reconstruction

   y=zeros(N,1);
   Lp=min(L,K);
   Kp=max(L,K);

   for k=0:Lp-2
     for m=1:k+1;
      y(k+1)=y(k+1)+(1/(k+1))*rca(m,k-m+2);
     end
   end

   for k=Lp-1:Kp-1
     for m=1:Lp;
      y(k+1)=y(k+1)+(1/(Lp))*rca(m,k-m+2);
     end
   end

   for k=Kp:N
      for m=k-Kp+2:N-Kp+1;
       y(k+1)=y(k+1)+(1/(N-k))*rca(m,k-m+2);
     end
   end

   figure;subplot(2,1,1);hold on;xlabel('Data poit');ylabel('Original and reconstructed series')
   plot(x1);grid on;plot(y,'r')

   r=x1-y;
   subplot(2,1,2);plot(r,'g');xlabel('Data poit');ylabel('Residual series');grid on
   vr=(sum(d(I))/sev)*100;


text:
% Chapter 17 - Neural Networks.
% Programs_17a - The Generalized Delta Learning Rule (Figure 17.7).
% Copyright Birkhauser 2013. Stephen Lynch.

function Programs_17a
% Load Boston housing data.
load housing.txt
X = housing(:,[6 9 13]);
t = housing(:,14);

% Scale to zero mean, unit variance and introduce bias on input.
xmean = mean(X);
xstd = std(X);
X = (X-ones(size(X,1),1)*xmean)./(ones(size(X,1),1)*xstd);
X = [ones(size(X,1),1) X];
tmean = (max(t)+min(t))/2;
tstd = (max(t)-min(t))/2;
t = (t-tmean)/tstd;

% Initialise random weight vector.
randn('seed', 123456);
w(:,1) = 0.1*randn(size(X,2),1);
y1 = tanh(X*w(:,1));
e1 = t-y1;
mse(1) = var(e1);

% Do numEpochs iterations.
numEpochs = 10;
numPatterns = size(X,1);
eta = 0.001;
k = 1;
for m=1:numEpochs
  for n=1:numPatterns
	% Calculate feedforward output, error, and gradient.
	yk = tanh(X(n,:)*w(:,k));
	err = yk-t(n);
	g = X(n,:)'*((1-yk^2)*err);
	% Update the weight.
	w(:,k+1) = w(:,k) - eta*g;
	k = k+1;
  end
end

for m=1:size(w,1)
plot(1:k, w(m,:))
hold on
end

fsize=15;
set(gca,'XTick',0:2000:6000,'FontSize',fsize)
set(gca,'YTick',-0.3:0.1:0.3,'FontSize',fsize)
xlabel('Number of Iterations','FontSize',fsize)
ylabel('Weights','FontSize',fsize)
hold off

% End of Programs_17a.


text:
% [p,inl,res,L] = ransacfit(x,model,maxres,smpls,smpln[,doLO,cnstr]) - Ransac fitting
%
% x         ... data - structure depends on the model, see the function
%                      implementation
% model     ... {'line','line-cnstr-seg','plane','conic','3Dpts','E5ptNister','P3PSolver','uX2PSolver'}
% maxres    ... maximal residual to vote
% smpls     ... sample size 0    = minimal sample
%                          (0,1) = percentile from the size of data
%                           >=1  = number of point in the samples
% smpln     ... number of samples
% doLO      ... {0,1} do local optimization (0 implicit)
% cnstr     ... constraint parameters
%               'line-cnstr-seg' ... cnstr = [x1 x2], such that x1 resp. x2
%                                                     is on the left resp. right
%                                                     from the line
% p         ... model parameters
% inl       ... inliers
% res       ... residuals for the model p
% L         ... fitting structure L{i} = output at the i-th improvement
%
% Example: >>ransacfit;

% T.Pajdla, pajdla@cmp.felk.cvut.cz
% 2015-07-11
function [p,inl,r,L] = at_ransac_p3p(x,maxres,smpls,smpln,doLO,bop)
if nargin>0 %% Fit
  if nargin<6, cnstr = []; end
  if nargin<5, doLO = false; end % implicitly do not use the local optimization
  %% model independent sample size
  % percentile from the data lengths
  if smpls >0 && smpls <1, smpl = max(2,ceil(length(x)*smpls)); end
  %         case 'P3PSolver'
  % x is struct containing
  % x.K = camera calibration matrix
  % x.iK = inversion of the camera calibration matris
  % x.x = [u;X] ... stacked hom coordis of image points & coords of 3D points
  if smpls == 0, smpl = 3; end % the minimal case
  if smpls>=3, smpl = smpls; end  % enforce the number of samples
  fitF = @P3PSolver; % 4pt absolute pose
  resF = @PerspRepErr; % residual function
  fitDataF = @P3PSolverFitImPoints; % original data
  resDataF = @E5ptResImPoints; % original data
  cnstr = x.K;

  % check if the sample size has been assigned
  if ~exist('smpl','var'), error('ransacfit: invalid sample size for smpls = %s',smpls); end
  %% RANSAC
  % prepare method depndent data from fitting and evaluation
  xf = fitDataF(x);
  xr = resDataF(x);
  % do ransac
  si = round((size(xf,2)-1)*rand(smpl,smpln)+1); % samples in the columns
  iN = 0;
  p = [];
  inl = [];
  r = [];
  k = 1;
  for i = si
    y = xf(:,i); % get a sample
    pf = fitF(y); % fit a model
    r  = resF(pf,xr,cnstr); % eval residuals & constraints
    il = abs(r)<=maxres; % inliers
    if size(il,1)>1 % there are more alternative models returned by fitF
      in = sum(il,2); % inlier #
      [in,ix] = max(in);
      il = il(ix,:); % select the best inliers
      r = r(ix,:); % select the residuals
      pf = pf{ix}; % select the best model
    else
      in = sum(il); % inlier #
      if ~isempty(pf)
        if iscell(pf)
          pf = pf{1};
        end
      end
    end
    if in>iN % larger support
      if doLO % local optimization
        y = xf(:,il); % get all inliers

%         po = fitF(y); % fit a model

        bop.constant_points = ones(1,sum(il));
        [uout,Kout,Rout,Cout,Xout,eout] = uPXBA_({y(1:2,:)},{pf},y(4:6,:),bop);
        po{1} = [Rout{1} -Rout{1}*Cout{1}];

        ro = resF(po,xr,cnstr); % eval residuals & constraints
        ilo = abs(ro)<=maxres; % inliers
        if size(ilo,1)>1 % there are more alternative models returned by fitF
          ino = sum(ilo,2); % inlier #
          [ino,ix] = max(ino);
          ilo = ilo(ix,:); % select the best inliers
          ro = ro(ix,:); % select the residuals
          po = po{ix}; % select the best model
        else
          ino = sum(ilo); % inlier #
        end
        if ino>in % an improvement
          p = po; r = ro; il = ilo; in = ino;
          fprintf('LO succeeded!\n');
        end
      end
      iN  = in;
      inl = il;
      p = pf;
      if nargout>3 % store the history
        L{k}.iN = iN; L{k}.p = p; L{k}.res = r; L{k}.inl = inl;
        k = k + 1;
      end
    end
  end
  % refit to all inliers (must return only one model!)
  if sum(inl)>0
    y   = xf(:,inl); % inliers
    p = fitF(y,p); % fit
    r = resF(p,xr,cnstr); % residuals
    inl = abs(r)<=maxres; % inliers
  else

  end
  if nargout>3 % store the history
    L{k}.iN = iN; L{k}.p = p; L{k}.res = r; L{k}.inl = inl;
  end

end

author: HajimeTaira


text:
%ComputeMarginal Computes the marginal over a set of given variables
%   M = ComputeMarginal(V, F, E) computes the marginal over variables V
%   in the distribution induced by the set of factors F, given evidence E
%
%   M is a factor containing the marginal over variables V
%   V is a vector containing the variables in the marginal e.g. [1 2 3] for
%     X_1, X_2 and X_3.
%   F is a vector of factors (struct array) containing the factors
%     defining the distribution
%   E is an N-by-2 matrix, each row being a variable/value pair.
%     Variables are in the first column and values are in the second column.
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function M = ComputeMarginal(V, F, E)

% Check for empty factor list
assert(numel(F) ~= 0, 'Error: empty factor list');

  F = ObserveEvidence(F, E);
  Joint = ComputeJointDistribution(F);
  Joint.val = Joint.val ./ sum(Joint.val);
  M = FactorMarginalization(Joint, setdiff(Joint.var, V));
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [gd,relres,iter] = gabconvexopt(g,a,M,varargin)
%GABCONVEXOPT Compute a window using convex optimization
%   Usage: gout=gabconvexopt(g,a,M);
%          gout=gabconvexopt(g,a,M, varagin);
%
%   Input parameters:
%     g      : Window function /initial point (tight case)
%     a      : Time shift
%     M      : Number of Channels
%
%   Output parameters:
%     gout   : Output window
%     iter   : Number of iterations
%     relres : Reconstruction error
%
%   `gabconvexopt(g,a,M)` computes a window *gout* which is the optimal
%   solution of the convex optimization problem below
%
%   .. gd  = argmin_x    || alpha x||_1 +  || beta Fx||_1
%
%   ..                 + || omega (x -g_l) ||_2^2 + delta || x ||_S0
%
%   ..                 + gamma || nabla F x ||_2^2 + mu || nabla x ||_2^2
%
%   ..     such that  x satifies the constraints
%
%   .. math:: \begin{split}  \text{gd}  = & \text{arg} \min_x    \| \alpha x \|_1 +  \| \beta \mathcal{F}x\|_1  \\ &  + \| \omega (x - g_l) \|_2^2  \\ &  \delta \| x \|_{S0}+ \mu \| \nabla x \|_2^2 +\gamma \| \nabla \mathcal{F} x \|_2^2 \\ &    \text{such that } x \text{ satisfies the constraints} \end{split}
%
%   Three constraints are possible:
%
%   * *x* is dual with respect of g
%
%   * *x* is tight
%
%   * *x* is compactly supported on Ldual
%
%   **Note**: This function require the unlocbox. You can download it at
%   `<http://unlocbox.sourceforge.net>`_
%
%   The function uses an iterative algorithm to compute the approximate.
%   The algorithm can be controlled by the following flags:
%
%     'alpha',alpha  Weight in time. If it is a scalar, it represent the
%                  weights of the entire L1 function in time. If it is a
%                  vector, it is the associated weight assotiated to each
%                  component of the L1 norm (length: Ldual).
%                  Default value is $\alpha=0$.
%                  **Warning**: this value should not be too big in order to
%                  avoid the the L1 norm proximal operator kill the signal.
%                  No L1-time constraint: $\alpha=0$
%
%     'beta',beta  Weight in frequency. If it is a scalar, it represent the
%                  weights of the entire L1 function in frequency. If it is a
%                  vector, it is the associated weight assotiated to each
%                  component of the L1 norm in frequency. (length: Ldual).
%                  Default value is $\beta=0$.
%                  **Warning**: this value should not be too big in order to
%                  avoid the the L1 norm proximal operator kill the signal.
%                  No L1-frequency constraint: $\beta=0$
%
%     'omega',omega  Weight in time of the L2-norm. If it is a scalar, it represent the
%                  weights of the entire L2 function in time. If it is a
%                  vector, it is the associated weight assotiated to each
%                  component of the L2 norm (length: Ldual).
%                  Default value is $\omega=0$.
%                  No L2-time constraint: $\omega=0$
%
%     'glike',g_l  $g_l$ is a windows in time. The algorithm try to shape
%                  the dual window like $g_l$. Normalization of $g_l$ is done
%                  automatically. To use option omega should be different
%                  from 0. By default $g_d=0$.
%
%     'mu', mu     Weight of the smooth constraint Default value is 1.
%                  No smooth constraint: $\mu=0$
%
%     'gamma', gamma  Weight of the smooth constraint in frequency. Default value is 1.
%                  No smooth constraint: $\gamma=0$
%
%     'delta', delta  Weight of the S0-norm. Default value is 0.
%                  No S0-norm: $\delta=0$
%
%     'support' Ldual  Add a constraint on the support. The windows should
%                  be compactly supported on Ldual.
%
%     'tight'      Look for a tight windows
%
%     'dual'       Look for a dual windows (default)
%
%     'painless'   Construct a starting guess using a painless-case
%                  approximation. This is the default
%
%     'zero'       Choose a starting guess of zero.
%
%     'rand'       Choose a random starting phase.
%
%     'tol',t      Stop if relative residual error is less than the
%                  specified tolerance.
%
%     'maxit',n    Do at most n iterations. default 200
%
%     'print'      Display the progress.
%
%     'debug'      Display all the progresses.
%
%     'quiet'      Don't print anything, this is the default.
%
%     'fast'       Fast algorithm, this is the default.
%
%     'slow'       Safer algorithm, you can try this if the fast algorithm
%                  is not working. Before using this, try to iterate more.
%
%     'printstep',p  If 'print' is specified, then print every p'th
%                    iteration. Default value is p=10;
%
%     'hardconstraint' Force the projection at the end (default)
%
%     'softconstaint' Do not force the projection at the end
%
%   See also: gaboptdual, gabdual, gabtight, gabfirtight, gabopttight

% Author: Nathanael Perraudin
% Date  : 18 Feb 2014

if nargin<4
  error('%s: Too few input parameters.',upper(mfilename));
end;

if numel(g)==1
  error('g must be a vector (you probably forgot to supply the window function as input parameter.)');
end;

definput.keyvals.L=[];
definput.keyvals.lt=[0 1];
definput.keyvals.tol=1e-6;
definput.keyvals.maxit=200;
definput.keyvals.printstep=10;
definput.flags.print={'quiet','print','debug'};
definput.flags.algo={'fast','slow'};
definput.flags.constraint={'hardconstraint','softconstaint'};
definput.flags.startphase={'painless','zero','rand'};
definput.flags.type={'dual','tight'};

definput.keyvals.alpha=0;
definput.keyvals.omega=0;
definput.keyvals.beta=0;
definput.keyvals.mu=1;
definput.keyvals.gamma=1;
definput.keyvals.vart=0;
definput.keyvals.varf=0;
definput.keyvals.var2t=0;
definput.keyvals.var2f=0;
definput.keyvals.support=0;
definput.keyvals.delta=0;
definput.keyvals.deltaw=0;
definput.keyvals.glike=zeros(size(g));

[flags,kv]=ltfatarghelper({'L','tol','maxit'},definput,varargin);

% Determine the window. The window /must/ be an FIR window, so it is
% perfectly legal to specify L=[] when calling gabwin
[g,info]=gabwin(g,a,M,[],kv.lt,'callfun',upper(mfilename));

if kv.support
    Ldual=kv.support;
    % Determine L. L must be longer than L+Ldual+1 to make sure that no convolutions are periodic
    L=dgtlength(info.gl+Ldual+1,a,M);
else
    L=length(g);
    Ldual=L;
end

b=L/M;

% Determine the initial guess
if flags.do_zero
  gd_initial=zeros(Ldual,1);
end;

if flags.do_rand
  gd_initial=rand(size(g));
end;

if flags.do_painless
  gsmall=long2fir(g,M);
  gdsmall=gabdual(gsmall,a,M);
  gd_initial=fir2long(gdsmall,Ldual);
end;

% -------- do the convex optimization stuff

% Define the long original window
glong=fir2long(g,L);

%gabframebounds(g,a,M)

% Initial point
xin=gd_initial;
xin=fir2long(xin,L);

% -- * Setting the different prox for ppxa *--
% ppxa will minimize all different proxes

% value test for the selection constraint
nb_priors=0;

% - variance -
    if kv.vart % constraint in time
        if flags.do_debug
            param_l1.verbose=1; % display the results
        else
            param_l1.verbose=0; % do not display anything
        end

         % alpha is a scalar
        if mod(L,2)
             w=[0:1:(L-1)/2,(L-1)/2:-1:1]';
        else
             w=[0:1:L/2-1,L/2:-1:1]';
        end
        w=w.^2/L;

        param_l1.weights=w;
        nb_priors=nb_priors+1;
        g11.prox= @(x,T) prox_l1(x,kv.vart*T,param_l1); % define the prox_l1 as operator
        g11.eval= @(x) kv.vart*norm(w.*x,1); % the objectiv function is the l1 norm
    else % no L1 in time constraint
        g11.prox= @(x,T) x;
        g11.eval= @(x) 0;
    end

% - variance -
    if kv.varf % constraint in time

        param_l1_fourier.A= @(x) 1/sqrt(L)*fft(x); % Fourier operator
        param_l1_fourier.At= @(x) sqrt(L)*ifft(x); % adjoint of the Fourier operator
        if flags.do_debug
            param_l1_fourier.verbose=1; % display the results
        else
            param_l1_fourier.verbose=0; % do not display anything
        end

        if mod(L,2)
             w=[0:1:(L-1)/2,(L-1)/2:-1:1]';
        else
             w=[0:1:L/2-1,L/2:-1:1]';
        end
        w=w.^2/L;

        param_l1_fourier.weights=w;
        nb_priors=nb_priors+1;
        g12.prox= @(x,T) prox_l1(x,kv.varf*T,param_l1_fourier); % define the prox_l1 as operator
        g12.eval= @(x) kv.varf*norm(w.*x,1); % the objectiv function is the l1 norm
    else % no L1 in time constraint
        g12.prox= @(x,T) x;
        g12.eval= @(x) 0;
    end

% - variance2 -
    if kv.var2t % constraint in time
        if flags.do_debug
            param_l2.verbose=1; % display the results
        else
            param_l2.verbose=0; % do not display anything
        end

         % alpha is a scalar
        if mod(L,2)
             w=[0:1:(L-1)/2,(L-1)/2:-1:1]';
        else
             w=[0:1:L/2-1,L/2:-1:1]';
        end
        w=w/sqrt(L);

        param_l2.weights=w;
        nb_priors=nb_priors+1;
        g13.prox= @(x,T) prox_l2(x,kv.var2t*T,param_l2); % define the prox_l1 as operator
        g13.eval= @(x) kv.var2t*norm(w.*x,2)^2; % the objectiv function is the l1 norm
    else % no L1 in time constraint
        g13.prox= @(x,T) x;
        g13.eval= @(x) 0;
    end

% - variance2 -
    if kv.var2f % constraint in time

        param_l2_fourier.A= @(x) 1/sqrt(L)*fft(x); % Fourier operator
        param_l2_fourier.At= @(x) sqrt(L)*ifft(x); % adjoint of the Fourier operator
        if flags.do_debug
            param_l2_fourier.verbose=1; % display the results
        else
            param_l2_fourier.verbose=0; % do not display anything
        end

        if mod(L,2)
             w=[0:1:(L-1)/2,(L-1)/2:-1:1]';
        else
             w=[0:1:L/2-1,L/2:-1:1]';
        end
        w=w/sqrt(L);

        param_l2_fourier.weights=w;
        nb_priors=nb_priors+1;
        g14.prox= @(x,T) prox_l2(x,kv.var2f*T,param_l2_fourier); % define the prox_l1 as operator
        g14.eval= @(x) kv.var2f*norm(w.*x,2)^2; % the objectiv function is the l1 norm
    else % no L1 in time constraint
        g14.prox= @(x,T) x;
        g14.eval= @(x) 0;
    end

% - small L1 norm in coefficient domain -
    if kv.alpha % constraint in time
        if flags.do_debug
            param_l1.verbose=1; % display the results
        else
            param_l1.verbose=0; % do not display anything
        end

        if length(kv.alpha)==1 % alpha is a scalar
            kv.alpha=ones(size(xin))*kv.alpha;
        end
        param_l1.weights=kv.alpha;
        nb_priors=nb_priors+1;
        g1.prox= @(x,T) prox_l1(x,T,param_l1); % define the prox_l1 as operator
        g1.eval= @(x) norm(kv.alpha.*x,1); % the objectiv function is the l1 norm
    else % no L1 in time constraint
        g1.prox= @(x,T) x;
        g1.eval= @(x) 0;
    end

% - small L1 norm in Fourier domain -
    if kv.beta %frequency constraint
        param_l1_fourier.A= @(x) 1/sqrt(L)*fft(x); % Fourier operator
        param_l1_fourier.At= @(x) sqrt(L)*ifft(x); % adjoint of the Fourier operator
        if flags.do_debug
            param_l1_fourier.verbose=1; % display the results
        else
            param_l1_fourier.verbose=0; % Do not display anything
        end

        if length(kv.beta)==1 % alpha is a scalar
            kv.beta=ones(size(xin))*kv.beta;
        end

        param_l1_fourier.weights=kv.beta;
        % Here are the step for the prox
        %   2) go into the Fourier domain (prox_l1)
        %   3) soft thresholding (prox_l1)
        %   4) back in the time domain (prox_l1)
        nb_priors=nb_priors+1;
        g3.prox= @(x,T) prox_l1(x,T,param_l1_fourier);

        g3.eval= @(x) norm(kv.beta.*fft(x),1); % objectiv function
    else % no L1 in frequency constraint
        g3.prox= @(x,T) x;
        g3.eval= @(x) 0; % objectiv function

    end

% - DUAL OR TIGHT?-
if flags.do_tight
    % tight windows
    g2.prox= @(x,T) gabtight(x,a,M); % set the prox
    g2.eval= @(x) norm(x-gabdual(x,a,M,L)); % objectiv function
else
% - projection on a B2 ball -
    % Frame-type matrix of the adjoint lattice
    %G=tfmat('dgt',glong,M,a);
    Fal=frame('dgt',glong,M,a);
    G=framematrix(Fal,L);
    d=[a/M;zeros(a*b-1,1)];

    % Using a B2 ball projection
    % || Gcut' x - b ||_2 < epsilon
%     param_proj.A = @(x) G'*x; % forward operator
%     param_proj.At = @(x) G*x; % adjoint operator
%     param_proj.y = d;
%     param_proj.maxit = 200;      % maximum of iteration
%     param_proj.tight=0;         % not a tight frame
%     param_proj.nu=norm(G)^2; % frame bound on Gcut'
%     param_proj.verbose=0;       % diplay summary at the end
%     param_proj.epsilon=10*eps;       % radius of the B2 ball
%     g2.prox= @(x,T) fast_proj_b2(x,T,param_proj); % set the prox

    % Using a direct projection (better solution)
    param_proj.verbose=flags.do_debug;
    param_proj.y=d;
    param_proj.A=G';
    param_proj.AAtinv=(G'*G)^(-1);
    g2.prox= @(x,T) proj_dual(x,T,param_proj); % set the prox
    g2.eval= @(x) norm(G'*x-d); % objectiv function
end

% SUPPORT CONSTRAINT
if kv.support
% - set null coefficient
    g4.prox = @(x,T) forceeven(fir2long(long2fir(x,Ldual),L));
    g4.eval = @(x) 0;

% - function apply the two projections thanks to a poc algorithm.
    if flags.do_tight

        G={g2,g4};
        paramPOCS.tol=20*eps;
        paramPOCS.maxit=5000;
        paramPOCS.verbose=flags.do_print+flags.do_debug;
        paramPOCS.abs_tol=1;
        g5.prox = @(x,T) pocs(x,G,paramPOCS);
        % g5.prox = @(x,T) ppxa(x,G,paramPOCS);
        % g5.prox = @(x,T) douglas_rachford(x,g2,g4,paramPOCS);
        % g5.prox = @(x,T) pocs2(x,g2,g4,20*eps,2000, flags.do_print+flags.do_debug);
        g5.eval = @(x) 0;

    else
        Fal=frame('dgt',glong,M,a);
        G=framematrix(Fal,L);
        d=[a/M;zeros(a*b-1,1)];
        Lfirst=ceil(Ldual/2);
        Llast=Ldual-Lfirst;
        Gcut=G([1:Lfirst,L-Llast+1:L],:);
        param_proj2.verbose=flags.do_debug;
        param_proj2.y=d;
        param_proj2.A=Gcut';
        param_proj2.AAtinv=pinv(Gcut'*Gcut);
        g5.prox= @(x,T) fir2long(proj_dual(long2fir(x,Ldual),T,param_proj2),L); % set the prox
        g5.eval= @(x) norm(G'*x-d); % objectiv function
    end

else
    g4.prox= @(x,T) x;
    g4.eval= @(x) 0; % objectiv function
    g5=g2;
end
% - function apply the two projections thanks to a douglas rachford algorithm.
%     param_douglas.verbose=1;
%     param_douglas.abs_tol=1;
%     param_douglas.maxit=2000;
%     param_douglas.tol=20*eps;
%     g6.prox = @(x,T) douglas_rachford(x,g2,g4,param_douglas);
%     g6.eval = @(x) 0;

% - small gradient norm -
% this is the smoothing parameter
    if kv.mu
        if flags.do_debug
            param_l2grad.verbose=1; % display the results
        else
            param_l2grad.verbose=0; % Do not display anything
        end
        nb_priors=nb_priors+1;
        g7.prox = @(x,T) prox_l2grad(fir2long(x,L),kv.mu*T,param_l2grad);
        g7.eval = @(x) norm(gradient(x))^2;
    else
        g7.prox = @(x,T) x;
        g7.eval = @(x) 0;
    end

% - small gradient norm in fourrier-
% this is the smoothing parameter
    if kv.gamma
        if flags.do_debug
            param_l2grad.verbose=1; % display the results
        else
            param_l2grad.verbose=0; % Do not display anything
        end
        nb_priors=nb_priors+1;
        g9.prox = @(x,T) prox_l2gradfourier(fir2long(x,L),kv.gamma*T,param_l2grad);
        g9.eval = @(x) norm(gradient(1/sqrt(L)*fft(x)))^2;
    else
        g9.prox = @(x,T) x;
        g9.eval = @(x) 0;
    end

% - small L2 norm in coefficient domain -
    if kv.omega % constraint in time
        if flags.do_debug
            param_l2.verbose=1; % display the results
        else
            param_l2.verbose=0; % do not display anything
        end

        if length(kv.omega)==1 % alpha is a scalar
            kv.alpha=ones(size(xin))*kv.omega;
        end
        param_l2.weights=kv.omega;
        if sum(kv.glike)
           kv.glike=fir2long(kv.glike,L);

           glike=kv.glike/norm(kv.glike)*norm(gabdual(g,a,M));
           param_l2.y=fir2long(glike,L);
        end
        nb_priors=nb_priors+1;
        g8.prox= @(x,T) prox_l2(x,T,param_l2); % define the prox_l2 as operator
        g8.eval= @(x) norm(kv.omega.*x-kv.glike,'fro')^2; % the objectiv function is the l2 norm
    else % no L1 in time constraint
        g8.prox= @(x,T) x;
        g8.eval= @(x) 0;
    end

% - small S0 norm -
    if kv.delta %frequency constraint
        gauss=pgauss(L,1);

        [A,B]=gabframebounds(gauss,1,L);
        AB=(A+B)/2;
        param_S0.A= @(x) dgt(x,gauss,1,L)/sqrt(AB);
        param_S0.At= @(x) idgt(x,gauss,1,L)/sqrt(AB);
        if flags.do_debug
            param_S0.verbose=1; % display the results
        else
            param_S0.verbose=0; % Do not display anything
        end

        nb_priors=nb_priors+1;
        g10.prox= @(x,T) prox_l1(x,T*kv.delta,param_S0);

        g10.eval= @(x) kv.delta*norm(reshape(dgt(x,gauss,1,L),[],1),1); % objectiv function
    else % no L1 in frequency constraint
        g10.prox= @(x,T) x;
        g10.eval= @(x) 0; % objectiv function

    end

% - small weighted S0 norm -
    if kv.deltaw %frequency constraint
        gauss=pgauss(L,1);

        [A,B]=gabframebounds(gauss,1,L);
        AB=(A+B)/2;
        param_S0.A= @(x) dgt(x,gauss,1,L)/sqrt(AB);
        param_S0.At= @(x) idgt(x,gauss,1,L)/sqrt(AB);
        if flags.do_debug
            param_S0.verbose=1; % display the results
        else
            param_S0.verbose=0; % Do not display anything
        end

        if mod(L,2)
             w=[0:1:(L-1)/2,(L-1)/2:-1:1]';
        else
             w=[0:1:L/2-1,L/2:-1:1]';
        end
        w=w/sqrt(L);

        %W=w*w';

        W=repmat(w,1,L).^2+repmat(w',L,1).^2;

        W=sqrt(W);

        param_S0.weights=W;
        nb_priors=nb_priors+1;
        g15.prox= @(x,T) prox_l1(x,T*kv.deltaw,param_S0);

        g15.eval= @(x) kv.deltaw*norm(reshape(dgt(x,gauss,1,L),[],1),1); % objectiv function
    else % no L1 in frequency constraint
        g15.prox= @(x,T) x;
        g15.eval= @(x) 0; % objectiv function

    end

% -- * PPXA function, the solver * --

% parameter for the solver
    param.maxit=kv.maxit; % maximum number of iteration
    param.tol=kv.tol;
    if flags.do_quiet
        param.verbose=0;
    end

    % Definition of the function f (the order is important)
    if flags.do_fast && flags.do_tight
        F={g1, g3,g7,g9,g8, g2, g4,g10,g11,g12,g13,g14,g15};

    else
        F={g1, g3,g7,g9,g8, g5,g10,g11,g12,g13,g14,g15};
    end

    % solving the problem

    if nb_priors
        [gd,iter,~]=ppxa(xin,F,param);

        % Force the hard constraint
        if flags.do_hardconstraint
            % In case of use of the douglas rachford algo instead of POCS
            %  gd=g6.prox(gd,0); % force the constraint

             gd=g5.prox(gd,0);
        end
    else
        fprintf( ' Warning!!! No prior selected! -- Only perform a projection. \n')
        gd=g5.prox(xin,0);
    end

    % compute the error
    if flags.do_tight
        relres=gabdualnorm(gd,gd,a,M,L);
    else
        relres=gabdualnorm(g,gd,a,M,L);
    end

   if kv.support
        % set the good size
        gd=long2fir(gd,Ldual);
   end

end

% function x=pocs2(x,g1,g2,tol,maxii,flagp)
% % this function implement a POCS algorithm, projection onto convex Set
% % using the differents projection of the algorithm.
% tola=1;
% ii=0;
% tola_old=tola;
% while (tola>tol)
%     x=g2.prox(g1.prox(x,0),0);
%     tola=g1.eval(x);
%     ii=ii+1;
%     if (logical(1-logical(mod(ii,50))) && flagp)
%        fprintf('      POCS sub-iteration: %i  -- Tol : %g\n',ii,tola)
%     end
%     if ii> maxii
%         break;
%     end
%     if abs(tola_old-tola)/tola<tol % avoid infinite loop
%         break;
%     end
%     tola_old=tola;
% end
% end

function x=forceeven(x)
% this function force the signal to be even
   x=  (x+involute(x))/2;
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% VL_DEMO  Run all demos

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

root = vl_root ;
demo_path = fullfile(root,'toolbox','demo') ;
test_path = fullfile(root,'toolbox','test') ;

addpath(demo_path) ;
addpath(test_path) ;

vl_demo_sift_basic ;
vl_demo_sift_peak ;
vl_demo_sift_edge ;
vl_demo_sift_cmd ;
vl_demo_sift_or ;
vl_demo_sift_match ;
vl_demo_sift_vs_ubc ;

vl_demo_dsift ;

vl_demo_mser_basic ;
vl_demo_mser_delta ;
vl_demo_mser_cmd ;

vl_test_hikmeans ;
vl_test_ikmeans ;
vl_test_imintegral ;

vl_demo_aib ;

vl_demo_quickshift ;
vl_demo_slic ;

vl_demo_kdtree ;
vl_demo_kdtree_sift ;
vl_demo_kdtree_self ;
vl_demo_kdtree_forest ;
vl_demo_kdtree_ann ;

vl_demo_imdisttf ;

rmpath(demo_path);
rmpath(test_path);

author: jianxiongxiao

repo: ProfXkit


text:
function [fevd_estimates]=fevd(struct_irf_record,gamma_record,It,Bu,n,IRFperiods,FEVDband)

% function [fevd_record]=fevd(struct_irf_record,gamma_record,It,Bu,IRFperiods,n)
% runs the gibbs sampler to obtain draws from the posterior distribution of FEVD
% inputs:  - cell 'struct_irf_record': record of the gibbs sampler draws for the orthogonalised IRFs
%          - matrix 'gamma_record': record of the gibbs sampler draws for the structural disturbances variance-covariance matrix gamma
%          - integer 'It': total number of iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'Bu': number of burn-in iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'IRFperiods': number of periods for IRFs
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
% outputs: - cell 'fevd_record': record of the gibbs sampler draws for the FEVD

% this function implements algorithm 3.1.1
% preliminary tasks
% define the time horizon of FEVD as that of IRFs
FEVDperiods=IRFperiods;

% create the first cell
temp=cell(n,n+1);

% now prepare the evaluation of (3.1.13)
% start by filling the first column of every Tij matrix in the cell
% loop over rows of temp
for jj=1:n
   % loop over columns of temp
   for ii=1:n
   % square each IRF element
   temp{jj,ii}(:,1)=struct_irf_record{jj,ii}(:,1).^2;
   end
end

% fill all the other entries of the Tij matrices
% loop over rows of temp
for jj=1:n
   % loop over columns of temp
   for ii=1:n
      % loop over remaining columns
      for kk=2:FEVDperiods
      % define the column as the square of the corresponding column in orthogonalised_irf_record
      % additioned to the value of the preceeding columns, which creates the cumulation
      temp{jj,ii}(:,kk)=struct_irf_record{jj,ii}(:,kk).^2+temp{jj,ii}(:,kk-1);
      end
   end
end

% multiply each matrix in the cell by the variance of the structural shocks
% to do so, loop over simulations (rows of the Tij matrices)
for jj=1:It-Bu
% recover the covariance matrix of structural shocks gamma for this iteration
gamma=reshape(gamma_record(:,jj),n,n);
% loop over rows of temp
   for ii=1:n
   % loop over columns of temp
      for kk=1:n
      % multiply row jj of the matrix by the variance of the structural shock
      temp{ii,kk}(jj,:)=temp{ii,kk}(jj,:)*gamma(kk,kk);
      end
   end
% then go for next iteration
end

% obtain now the values for Ti, the (n+1)th matrix of each row
% loop over rows of temp
for jj=1:n
% start the summation over Tij matrices
temp{jj,n+1}=temp{jj,1};
   % sum over remaining columns
   for ii=2:n
   temp{jj,n+1}=temp{jj,n+1}+temp{jj,ii};
   end
end

% create the output cell fevd_record
fevd_record=cell(n,n);

% fill the cell
% loop over rows of fevd_record
for jj=1:n
   % loop over columns of fevd_record
   for ii=1:n
   % define the matrix Vfij as the division (pairwise entry) of Tfij by Tfj
   fevd_record{jj,ii}=temp{jj,ii}./temp{jj,n+1};
   end
end

%% create the FEVD estimates output
% create first the cell that will contain the estimates
fevd_estimates=cell(n,n);

% for each variable and each variable contribution along with each period, compute the median, lower and upper bound from the Gibbs sampler records
% consider variables in turn
for ii=1:n
   % consider contributions in turn
   for jj=1:n
      % consider periods in turn
      for kk=1:IRFperiods
      % compute first the lower bound
      fevd_estimates{ii,jj}(1,kk)=quantile(fevd_record{ii,jj}(:,kk),(1-FEVDband)/2);
      % then compute the median
      fevd_estimates{ii,jj}(2,kk)=quantile(fevd_record{ii,jj}(:,kk),0.5);
      % finally compute the upper bound
      fevd_estimates{ii,jj}(3,kk)=quantile(fevd_record{ii,jj}(:,kk),1-(1-FEVDband)/2);
      end
   end
end

author: european-central-bank

repo: BEAR-toolbox