text:
function dlap_file_print ( n, nelt, isym, irhs, isoln, ia, ja, a, rhs, ...
  soln, title )

%*****************************************************************************80
%
%% DLAP_FILE_PRINT prints a DLAP linear system that was stored in a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Mark Seager,
%    A SLAP for the Masses,
%    Lawrence Livermore National Laboratory,
%    Technical Report UCRL-100267, December 1988.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer NELT, the number of non-zeros stored in A.
%
%    Input, integer ISYM, a flag to indicate symmetric
%    storage format.
%    * 0, all nonzero entries of the matrix are stored.
%    * 1, the matrix is symmetric, and only the lower triangle of the
%    matrix is stored.
%
%    Input, integer IRHS, is 1 if a right hand side vector
%    is included.
%
%    Input, integer ISOLN, is 1 if a solution vector is included.
%
%    Input, integer IA(NELT), integer JA(NELT),
%    real A(NELT), the DLAP triad matrix description.
%
%    Input, real RHS(N), the right hand side vector.
%
%    Input, real SOLN(N), the solution to the linear system.
%
%    Input, string TITLE, a title to be printed.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  dlap_header_print ( n, nelt, isym, irhs, isoln );
%
%  Write out the matrix.
%
  r8sp_print ( n, n, nelt, isym, ia, ja, a, '  The sparse matrix' );
%
%  Write the right hand side.
%
  if ( irhs == 1 )
    dlap_rhs_print ( n, rhs );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  No right hand side vector was supplied.\n' );
  end
%
%  Write the solution.
%
  if ( isoln == 1 )
    dlap_soln_print ( n, soln );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  No solution vector was supplied.\n' );
  end

  return
end


text:
classdef LvdPluginOptimVarSet < matlab.mixin.SetGet
    %LvdPluginOptimVarSet Summary of this class goes here
    %   Detailed explanation goes here

    properties
        pluginVars(1,:) LvdPluginOptimVarWrapper

        lvdData LvdData
    end

    methods
        function obj = LvdPluginOptimVarSet(lvdData)
            obj.lvdData = lvdData;
        end

        function addPluginVar(obj, newPluginVar)
            obj.pluginVars(end+1) = newPluginVar;
        end

        function removePluginVar(obj, pluginVar)
            obj.pluginVars([obj.pluginVars] == pluginVar) = [];
        end

        function [listBoxStr, pluginVars] = getListboxStr(obj)
            listBoxStr = {};

            for(i=1:length(obj.pluginVars)) %#ok<*NO4LP>
                if(obj.pluginVars(i).isVariableActive())
                    varStr = '* ';
                else
                    varStr = '';
                end

                listBoxStr{end+1} = sprintf('%u - %s%s', i, varStr, obj.pluginVars(i).name); %#ok<AGROW>
            end

            pluginVars = obj.pluginVars;
        end

        function pluginVars = getPluginVarsArray(obj)
            pluginVars = obj.pluginVars;
        end

        function inds = getIndsForPluginVars(obj, pluginVars)
            inds = find(ismember(obj.pluginVars, pluginVars));
        end

        function plugin = getPluginVarAtInd(obj, ind)
            if(ind > 0 && ind <= length(obj.pluginVars))
                plugin = obj.pluginVars(ind);
            else
                plugin = LvdPluginOptimVarWrapper.empty(1,0);
            end
        end

        function numPluginVars = getNumPluginVars(obj)
            numPluginVars = length(obj.pluginVar);
        end

        function varValues = getPluginVarValues(obj)
            varValues = [obj.pluginVars.value];
            varValues = varValues(:);
        end

        function tf = isVarAPluginVar(obj, var)
            if(not(isempty(obj)))
                if(isempty(obj.pluginVars))
                    tf = false;

                else
                    vars = [obj.pluginVars.optVar];
                    tf = any(var == vars);
                end

            else
                tf = true;
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function result=robpca(x,varargin)

%ROBPCA is a 'ROBust method for Principal Components Analysis'.
% It is resistant to outliers in the data. The robust loadings are computed
% using projection-pursuit techniques and the MCD method.
% Therefore ROBPCA can be applied to both low and high-dimensional data sets.
% In low dimensions, the MCD method is applied (see mcdcov.m).
%
% The ROBPCA method is described in
%  Hubert, M., Rousseeuw, P.J., Vanden Branden, K. (2005), ROBPCA: a
%  new approach to robust principal components analysis, Technometrics, 47, 64-79.
%
%
% To select the number of principal components, a robust PRESS (predicted
% residual sum of squares) curve is drawn, based on a fast algorithm for
% cross-validation. This approach is described in:
%
%    Hubert, M., Engelen, S. (2007),
%    "Fast cross-validation of high-breakdown resampling algorithms for PCA",
%    Computational Statistics and Data Analysis, 51, 5013-5024.
%
% ROBPCA is designed for normally distributed data. If the data are skewed,
% a modification of ROBPCA is available, based on the adjusted outlyingness
% (see adjustedoutlyingness.m). This method is described in:
%
%    Hubert, M., Rousseeuw, P.J., Verdonck, T. (2008),
%    "Robust PCA for skewed data and its outlier map", Computational Statistics
%    and Data Analysis, in press.
%
%  For the up-to-date reference, please consult the website:
%    wis.kuleuven.be/stat/robust.html
%
% Required input arguments:
%            x : Data matrix (observations in the rows, variables in the
%                columns)
%
% Optional input arguments:
%            k : Number of principal components to compute. If k is missing,
%                or k = 0, a scree plot and a press curve are drawn which allows you to select
%                the number of principal components.
%         kmax : Maximal number of principal components to compute (default = 10).
%                If k is provided, kmax does not need to be specified, unless k is larger
%                than 10.
%        alpha : (1-alpha) measures the fraction of outliers the algorithm should
%                resist. Any value between 0.5 and 1 may be specified (default = 0.75).
%            h : (n-h+1) measures the number of outliers the algorithm should
%                resist. Any value between n/2 and n may be specified. (default = 0.75*n)
%                Alpha and h may not both be specified.
%          mcd : If equal to one: when the number of variables is sufficiently small,
%                the loadings are computed as the eigenvectors of the MCD covariance matrix,
%                hence the function 'mcdcov.m' is automatically called. The number of
%                principal components is then taken as k = rank(x). (default)
%                If equal to zero, the robpca algorithm is always applied.
%        plots : If equal to one, a scree plot, a press curve and a robust score outlier map are
%                drawn (default). If the input argument 'classic' is equal to one,
%                the classical plots are drawn as well.
%                If 'plots' is equal to zero, all plots are suppressed (unless k is missing,
%                then the scree plot and press curve are still drawn).
%                See also makeplot.m
%        labsd : The 'labsd' observations with largest score distance are
%                labeled on the outlier map. (default = 3)
%        labod : The 'labod' observations with largest orthogonal distance are
%                labeled on the outlier map. default = 3)
%      classic : If equal to one, the classical PCA analysis will be performed
%                (see also cpca.m). (default = 0)
%        scree : If equal to one, a scree plot is drawn. If k is given as input, the default value is 0, else the default value is one.
%        press : If equal to one, a plot of robust press-values is drawn.
%                If k is given as input, the default value is 0, else the default value is one.
%                If the input argument 'skew' is equal to one, no plot is
%                drawn.
%    robpcamcd : If equal to one (default), the whole robpca procedure is run (computation of outlyingness and
%                MCD).
%                If equal to zero, the program stops after the computation of the outlyingness. The
%                robust eigenvectors then correspond with the eigenvectors of the covariance matrix
%                of the h observations with smallest outlyingness. This yields the same
%                PCA subspace as the full robpca, but not the same eigenvectors and eigenvalues.
%         skew : If equal to zero the regular robpca is run. If equal to
%                one, the adjusted robpca algorithm for skewed data is run.
%
% I/O: result=robpca(x,'k',k,'kmax',10,'alpha',0.75,'h',h,'mcd',1,'plots',1,'labsd',3,'labod',3,'classic',0);
%  The user should only give the input arguments that have to change their default value.
%  The name of the input arguments needs to be followed by their value.
%  The order of the input arguments is of no importance.
%
% Examples:
%    result=robpca(x,'k',3,'alpha',0.65,'plots',0)
%    result=robpca(x,'alpha',0.80,'kmax',15,'labsd',5)
%
% The output of ROBPCA is a structure containing
%
%    result.P        : Robust loadings (eigenvectors)
%    result.L        : Robust eigenvalues
%    result.M        : Robust center of the data
%    result.T        : Robust scores
%    result.k        : Number of (chosen) principal components
%    result.kmax     : Maximal number of principal components
%    result.alpha    : see interpretation in the list of input arguments
%    result.h        : The quantile h used throughout the algorithm
%    result.Hsubsets : A structure that contains H0, H1 and Hfreq:
%                      H0 : The h-subset that contains the h points with the smallest outlyingness.
%                      H1 : The optimal h-subset of mcdcov.
%                      Hfreq : The subset of h points which are the most frequently selected during the mcdcov
%                              algorithm.
%    result.sd       : Robust score distances within the robust PCA subspace
%    result.od       : Orthogonal distances to the robust PCA subspace
%    result.cutoff   : Cutoff values for the robust score and orthogonal distances
%    result.flag     : The observations whose score distance is larger than result.cutoff.sd (==> result.flag.sd)
%                      or whose orthogonal distance is larger than result.cutoff.od (==> result.flag.od)
%                      can be considered as outliers and receive a flag equal to zero (result.flag.all).
%                      The regular observations receive a flag 1.
%    result.class    : 'ROBPCA'
%    result.classic  : If the input argument 'classic' is equal to one, this structure
%                     contains results of the classical PCA analysis (see also cpca.m).
%
% Short description of the method:
%
% Let n denote the number of observations, and p the number of original variables,
% then ROBPCA finds a robust center (p x 1) of the data M and a loading matrix P which
% is (p x k) dimensional. Its columns are orthogonal and define a new coordinate
% system. The scores (n x k) are the coordinates of the centered observations with
% respect to the loadings: T=(X-M)*P.
% Note that ROBPCA also yields a robust covariance matrix (often singular) which
% can be computed as
%                         cov=out.P*out.L*out.P'
%
% To select the number of principal components, it is useful to look at the scree plot which
% shows the eigenvalues, and the press curve which displays a weighted sum of the squared
% cross-validated orthogonal distances.
% The outlier map visualizes the observations by plotting their orthogonal
% distance to the robust PCA subspace versus their robust distances
% within the PCA subspace. This allows to classify the data points into 4 types:
% regular observations, good leverage points, bad leverage points and
% orthogonal outliers.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust
%
% Written by Mia Hubert, Sabine Verboven, Karlien Vanden Branden, Sanne Engelen, Tim Verdonck
% Last Update: 17/06/2003,  03/07/2006,  31/07/2007
% Last Revision: 27/03/2008, 09/06/2008

%
% initialization with defaults
%
data=x;
[n,p]=size(data);

% First Step: classical PCA on data
if n < p
    [P1,T1,L1,r,Xc,clm]=kernelEVD(data);
else
	[P1,T1,L1,r,Xc,clm]=classSVD(data);
end
% dim(P1): p x r

if r==0
    error('All data points collapse!')
end

niter=100;
counter=1;
kmax=min([10,floor(n/2),r]);
k=0;
alfa=0.75;
h=min(floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*alfa),n);
labsd=3;
labod=3;
plots=1;
scree = 1;
press = 1;
mcd=1; % user wants the mcd approach (in case n>>p)
robpcamcd = 1;
cutoff  = 0.975;
skew=0;
% default is a structure needed for input checking
default=struct('alpha',alfa,'h',h,'labsd',labsd,'labod',labod,...
    'k',k,'plots',plots,'kmax',kmax,'mcd',mcd,'classic',0,'scree',scree,'press',press,'robpcamcd',robpcamcd,'cutoff',cutoff,'skew',0);
list=fieldnames(default);
options=default; %input by user
IN=length(list);
i=1;
%
if nargin==2
    error('Incorrect number of input arguments!')
end
dummy = 0;  %Assume we didn't get h or alpha, unless we find it below.
if nargin>2
    %
    %placing inputfields in array of strings
    %
    for j=1:nargin-2
        if rem(j,2)~=0
            chklist{i}=varargin{j};
            i=i+1;
        end
    end
    dummy=sum(strcmp(chklist,'h')+2*strcmp(chklist,'alpha')); %checking if h and alpha are provided both or not
    switch dummy
    case 0 % Take on default values
        options.alpha=alfa;
        if any(strcmp(chklist,'kmax'))
            for j=1:nargin-2 % searching the index of the accompanying field
                if rem(j,2)~=0 % fieldnames are placed on odd index
                    if strcmp('kmax',varargin{j})
                        I=j;
                    end
                end
            end
            options=setfield(options,'kmax',varargin{I+1});
            kmax = max(min([floor(options.kmax),floor(n/2),r]),1); %acceptable kmax
            options.h=min(floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*alfa),n); %depends on kmax, so if kmax is given by user, get it first!
        else %kmax is not given by user
           options.h=h;
        end
    case 3
        error('Both inputarguments alpha and h are provided. Only one is required.')
    end

    %
    % Checking which default parameters have to be changed
    % and keep them in the structure 'options'.
    %
    while counter<=IN
        index=strmatch(list(counter,:),chklist,'exact');
        if ~isempty(index) % in case of similarity
            for j=1:nargin-2 % searching the index of the accompanying field
                if rem(j,2)~=0 % fieldnames are placed on odd index
                    if strcmp(chklist{index},varargin{j})
                        I=j;
                    end
                end
            end
            options=setfield(options,chklist{index},varargin{I+1});
            index=[];
        end
        counter=counter+1;
    end
    options.h=floor(options.h);
    options.kmax=floor(options.kmax);
    options.k=floor(options.k);
    kmax=max(min([options.kmax,floor(n/2),r]),1);
    labod=max(0,min(floor(options.labod),n));
    labsd=max(0,min(floor(options.labsd),n));
    k=options.k;

    if k<0
        k=0;
    elseif k > kmax
        k=kmax;
        mess=sprintf(['Attention (robpca.m): The number of principal components, k = ',num2str(options.k)...
            ,'\n is larger than kmax= ',num2str(kmax),'; k is set to ',num2str(kmax)]);
        disp(mess)
    end
    if dummy==1 % checking input variable h
        options.alpha=options.h/n;
        if k==0
            if options.h < floor((n+kmax+1)/2 )
                options.h=floor((n+kmax+1)/2);
                options.alpha=options.h/n;
                mess=sprintf(['Attention (robpca.m): h should be larger than (n+kmax+1)/2.\n',...
                        'It is set to its minimum value ',num2str(options.h)]);
                disp(mess)
            end
        else
            if options.h < floor((n+k+1)/2)
                options.h=floor((n+k+1)/2);
                options.alpha=options.h/n;
                mess=sprintf(['Attention (robpca.m): h should be larger than (n+k+1)/2.\n',...
                        'It is set to its minimum value ',num2str(options.h)]);
                disp(mess)
            end
        end
        if options.h > n
            options.alpha=0.75;
            if k==0
                options.h=floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*options.alpha);
            else
                options.h=floor(2*floor((n+k+1)/2)-n+2*(n-floor((n+k+1)/2))*options.alpha);
            end
            mess=sprintf(['Attention (robpca.m): h should be smaller than n. \n',...
                    'It is set to its default value ',num2str(options.h)]);
            disp(mess)
        end
    elseif dummy==2 %checking input variable alpha
        if options.alpha < 0.5
            options.alpha=0.5;
            mess=sprintf(['Attention (robpca.m): Alpha should be larger than 0.5.\n',...
                    'It is set to 0.5.']);
            disp(mess)
        end
        if options.alpha > 1
            options.alpha=0.75;
            mess=sprintf(['Attention (robpca.m): Alpha should be smaller than 1. \n',...
                    'It is set to 0.75.']);
            disp(mess)
        end
        if k==0
            options.h=floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*options.alpha);
        else
            options.h=floor(2*floor((n+k+1)/2)-n+2*(n-floor((n+k+1)/2))*options.alpha);
        end
    end
    alfa=options.alpha;
    dummyh = strcmp(chklist,'h');
    dummykmax = strcmp(chklist,'kmax');
    %     if all(dummyh == 0) & any(dummykmax) & k==0
    %         h = min(floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*alfa),n);
    %     end
    if all(dummyh == 0)&& any(dummykmax) %kmax was given by the user
        if k==0
            options.h=floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*options.alpha);
        else
            options.h=floor(2*floor((n+k+1)/2)-n+2*(n-floor((n+k+1)/2))*options.alpha);
        end
    elseif all(dummyh == 0) && ~any(dummykmax) %kmax is the default value
        if k==0
            options.h=floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*options.alpha);
        else
            options.h=floor(2*floor((n+k+1)/2)-n+2*(n-floor((n+k+1)/2))*options.alpha);
        end
    end
    h=options.h;
    dummyscree = strcmp(chklist,'scree');
    dummypress = strcmp(chklist,'press');
    if all(dummyscree == 0)
        if k~=0
            options.scree = 0;
        end
    end
    if all(dummypress == 0)
        if k~=0
            options.press = 0;
        end
    end
    scree = options.scree;
    press = options.press;
    labsd=floor(max(0,min(options.labsd,n)));
    labod=floor(max(0,min(options.labod,n)));
    plots=options.plots;
    mcd=options.mcd;
    robpcamcd = options.robpcamcd;
    cutoff = options.cutoff;
    skew = options.skew;
    if skew==1
        press=0; %press curve not available for skewed data
    end
end
%
% MAIN PART
%
X=T1;
center=clm;
rot=P1;
% Depending on n and p, perform MCD or ROBPCA:
% p << n => MCD
p1=size(X,2);
if p1<=min(floor(n/5),kmax) && mcd && (skew==0)
    options.h=h;
    [res,raw]=mcdcov(X,'h',h,'plots',0);
    [U,S,P]=svd(res.cov,0);
    L=diag(S);
    if k~=0
        options.k=min(k,p1);
    else
        bdwidth=5;
        topbdwidth=30;
        set(0,'Units','pixels');
        scnsize=get(0,'ScreenSize');
        pos1=[bdwidth, 1/3*scnsize(4)+bdwidth, scnsize(3)/2-2*bdwidth, scnsize(4)/2-(topbdwidth+bdwidth)];
        pos2=[pos1(1)+scnsize(3)/2, pos1(2), pos1(3), pos1(4)];
        if press == 1
            outcvMcd = cvMcd(X,p1,res,h);
            figure('Position',pos1)
            set(gcf,'Name', 'PRESS curve','NumberTitle', 'off');
            plot(1:p1,outcvMcd.press,'o-')
            title('MCD')
            xlabel('number of LV')
            ylabel('R-PRESS')
        end
        if scree == 1
            figure('Position',pos2)
            screeplot(L,'MCD');
        end
        if (scree == 1) || (press == 1)
            cumperc = cumsum(L)./sum(L);
            disp(['The cumulative percentage of variance explained by the first ',num2str(kmax),' components is:']);
            disp(num2str(cumperc'));
            disp(['How many principal components would you like to retain? Max = ',num2str(kmax),'. ']);
            k=input('');
        end
        % to close the figures.
        if scree == 1
            close
        end
        if press == 1
            close
        end
    end
    options.k = k;
    T=(X-repmat(res.center,size(X,1),1))*U;
    out.M=center+res.center*rot';
    out.L=L(1:options.k)';
    out.P=rot*U(:,1:options.k);
    out.T=T(:,1:options.k);
    out.h=h;
    out.k=options.k;
    out.alpha=alfa;
    out.Hsubsets.H0 = res.Hsubsets.Hopt;
    out.Hsubsets.H1 = [];
    out.Hsubsets.Hfreq = res.Hsubsets.Hfreq;
    out.skew=skew;
else
    % p > n => ROBPCA
    niter=100;
    seed=0;
    if h~=n
        if skew==0
            B=twopoints(T1,250,seed); %n*ri
            for i=1:size(B,1)
                Bnorm(i)=norm(B(i,:),2);
            end
            Bnormr=Bnorm(Bnorm > 1.e-12); %ndirect*1
            B=B(Bnorm > 1.e-12,:);       %ndirect*n
            A=diag(1./Bnormr)*B;         %ndirect*n
            %projected points in columns
            Y=T1*A';%n*ndirect
            m=length(Bnormr);
            Z=zeros(n,m);
            for i=1:m
                [tmcdi,smcdi,weights]=unimcd(Y(:,i),h);
                if smcdi<1.e-12
                    r2=rank(data(weights,:));
                    if r2==1
                        error(['At least ',num2str(sum(weights)),' obervations are identical.']);
                    end
                else
                    Z(:,i)=abs(Y(:,i)-tmcdi)/smcdi;
                end
            end
            d=max(Z,[],2);
        else %adjusted robpca for skewed data
            outAO=adjustedoutlyingness(T1,'ndir',min(250*p,2500));
            d=outAO.adjout;
        end
        [ds,is]=sort(d);
        Xh=T1(is(1:h),:); % Xh contains h (good) points out of Xcentr
        [P2,T2,L2,r2,Xm,clmX]=classSVD(Xh);
        out.Hsubsets.H0 = is(1:h);
        Tn=(T1-repmat(clmX,n,1))*P2;
    else
        P2=eye(r);
        Tn=T1;
        L2=L1;
        r2=r;
        out.Hsubsets.H0=1:n;
        Xm=T1;
        clmX=zeros(1,size(T1,2));
    end

    %dim(P2) = r x r2
    L=L2;
    kmax=min(r2,kmax);

    % choice of k:
    %-------------
    bdwidth=5;
    topbdwidth=30;
    set(0,'Units','pixels');
    scnsize=get(0,'ScreenSize');
    pos1=[bdwidth, 1/3*scnsize(4)+bdwidth, scnsize(3)/2-2*bdwidth, scnsize(4)/2-(topbdwidth+bdwidth)];
    pos2=[pos1(1)+scnsize(3)/2, pos1(2), pos1(3), pos1(4)];

    if press == 1
        disp('The robust press curve based on cross-validation is now computed.')
        outprMCDkmax = projectMCD(Tn,L,kmax,h,niter,rot,P1,P2,center,cutoff);
        if size(out.Hsubsets.H0,2)==1
            out.Hsubsets.H0=out.Hsubsets.H0';
        end
        outprMCDkmax.Hsubsets.H0 = out.Hsubsets.H0;
        outpress = cvRobpca(data,kmax,outprMCDkmax,0,h);
        figure('Position',pos1)
        set(gcf,'Name', 'PRESS curve','NumberTitle', 'off');
        plot(1:kmax,outpress.press,'o-')
        title('ROBPCA')
        xlabel('Number of LV')
        ylabel('R-PRESS')
    else
        if size(out.Hsubsets.H0,2)==1
            out.Hsubsets.H0=out.Hsubsets.H0';
        end
    end

    if scree == 1
        figure('Position',pos2)
        screeplot(L(1:kmax),'ROBPCA')
    end

    if (scree == 1)||(press == 1)
        cumperc = (cumsum(L(1:kmax))./sum(L))';
        disp(['The cumulative percentage of variance explained by the first ',num2str(kmax),' components is:']);
        disp(num2str(cumperc));
        disp(['How many principal components would you like to retain? Max = ',num2str(kmax),'.']);
        k=input('');
        k=max(min(min(r2,k),kmax),1);
        % we compute again the robpca results for a specific k value. alpha
        % and h can change again, because until now they were based on the kmax
        % value.
        if dummy == 2
            options.h=floor(2*floor((n+k+1)/2)-n+2*(n-floor((n+k+1)/2))*alfa);
            %if dummy == 1 no changes needed
        elseif dummy~=1
            options.h=floor(2*floor((n+k+1)/2)-n+2*(n-floor((n+k+1)/2))*options.alpha);
        end
        h=options.h;
        % to close the figures.
        if scree == 1
            close
        end
        if press == 1
            close
        end
    else
        k=min(min(r2,k),kmax);
    end

    if k~=r % extra reweighting step
        XRc=T1-repmat(clmX,n,1);
        Xtilde=XRc*P2(:,1:k)*P2(:,1:k)';
        Rdiff = XRc-Xtilde;
        for i=1:n
            odh(i,1)=norm(Rdiff(i,:));
        end
        if skew==0
            [m,s]=unimcd(odh.^(2/3),h);
            cutoffodh = sqrt(norminv(cutoff,m,s).^3);
        else %adjusted robpca for skewed data
            mcodh=mc(odh);
            if mcodh>0
                cutoffodh = prctile(odh,75)+1.5*exp(3*mcodh)*iqr(odh);
            else
                cutoffodh = prctile(odh,75)+1.5*iqr(odh);
            end
            ttup = sort(-odh(odh<cutoffodh));
            cutoffodh = -ttup(1);
        end
        indexset = find(odh<=cutoffodh)';
        [P2,Th,Lh,rh,Xm,clmX]=classSVD(T1(indexset,:));
        if k>rh
            k = rh;
        end
    end
    center=center+clmX*rot';
    rot=rot*P2(:,1:k);
    Tn=(T1-repmat(clmX,n,1))*P2;

    % if only the subspace is important, not the PC themselves: do not
    % perform MCD anymore.
    if ~robpcamcd
        out.P = rot; %=P1*P2(:,1:k);
        out.T = Tn(:,1:k);
        out.M = center;
        out.L = Lh;
        out.k = k;
        out.h = h;
        out.alpha = alfa;
        out.kmax=kmax;
        out.skew=skew;
    end

    % projection, mcd
    %-----------------
    if skew==0
        outpr = projectMCD(Tn,L,k,h,niter,rot,P1,P2,center,cutoff);
    else %adjusted robpca for skewed data
        outpr = projectAO(Tn,k,h,rot,center);
    end
    out.T = outpr.T;
    out.P = outpr.P;
    out.M = outpr.M;
    out.L = outpr.L;
    out.k = k;
    out.kmax=kmax;
    out.h = h;
    out.alpha = alfa;
    out.skew = skew;
    if skew==0
        out.Hsubsets.H1 = outpr.Hsubsets.H1;
        out.Hsubsets.Hfreq = outpr.Hsubsets.Hfreq;
    else
        out.AO=outpr.AO;
        out.cutoff.AO=outpr.cutoff.AO;
    end
end

% Classical analysis
if options.classic==1
    out.classic.P=P1(:,1:out.k);
    out.classic.L=L1(1:out.k)';
    out.classic.M=clm;
    out.classic.T=T1(:,1:out.k);
    out.classic.k=out.k;
    out.classic.Xc=Xc;
end

outpr = out;

% Calculation of the distances, flags
%-------------------------------------

if options.classic == 1
    outDist = CompDist(data,r,outpr,cutoff,robpcamcd,out.classic);
else
    outDist = CompDist(data,r,outpr,cutoff,robpcamcd);
end

out.sd = outDist.sd;
out.cutoff.sd = outDist.cutoff.sd;
out.od = outDist.od;
out.cutoff.od = outDist.cutoff.od;
out.flag = outDist.flag;
out.class = outDist.class;
out.classic = outDist.classic;
if options.classic == 1
    out.classic.sd = outDist.classic.sd;
    out.classic.od = outDist.classic.od;
    out.classic.cutoff.sd = outDist.classic.cutoff.sd;
    out.classic.cutoff.od = outDist.classic.cutoff.od;
    out.classic.class = outDist.classic.class;
    out.classic.flag = outDist.classic.flag;
end

result=struct('P',{out.P},'L',{out.L},'M',{out.M},'T',{out.T},'k',{out.k},'kmax',{kmax},'alpha',{out.alpha},...
    'h',{out.h},'Hsubsets',{out.Hsubsets},'sd', {out.sd},'od',{out.od},'cutoff',{out.cutoff},'flag',out.flag',...
    'class',{out.class},'classic',{out.classic});

% Plots
try
    if plots && options.classic
        makeplot(result,'classic',1,'labsd',labsd,'labod',labod)
    elseif plots
        makeplot(result,'labsd',labsd,'labod',labod)
    end
catch %output must be given even if plots are interrupted
    %> delete(gcf) to get rid of the menu
end

%--------------------------------------------------------------------------
function outprMCD = projectMCD(Tn,L,k,h,niter,rot,P1,P2,center,cutoff)

% this function performs the last part of ROBPCA when k is determined.
% input :
%   Tn : the projected data
%   L  : the matrix of the eigenvalues
%   k  : the number of components
%   h  : lower bound for regular observations
%   niter : the number of iterations
%   rot : the rotation matrix
%   P1, P2: the different eigenvector matrices after each transformation
%   center : the classical center of the data

X2=Tn(:,1:k);
n = size(X2,1);
rot=rot(:,1:k);
% first apply c-step with h points from first step,i.e. those that
% determine the covariance matrix after the c-steps have converged.
mah=libra_mahalanobis(X2,zeros(size(X2,2),1),'cov',L(1:k));
oldobj=prod(L(1:k));
P4=eye(k);
korig=k;
for j=1:niter
    [mahs,is]=sort(mah);
    Xh=X2(is(1:h),:);
    [P,T,L,r3,Xm,clmX]=classSVD(Xh);
    obj=prod(L);
    X2=(X2-repmat(clmX,n,1))*P;
    center=center+clmX*rot';
    rot=rot*P;
    mah=libra_mahalanobis(X2,zeros(size(X2,2),1),'cov',diag(L));
    P4=P4*P;
    if ((r3==k) && (abs(oldobj-obj) < 1.e-12))
        break;
    else
        oldobj=obj;
        j=j+1;
        if r3 < k
            j=1;
            k=r3;
        end
    end
end
% dim(P4): k x k0 with k0 <= k but denoted as k
% dim X2: n x k0
% perform mcdcov on X2
[zres,zraw]= mcdcov(X2,'plots',0,'ntrial',250,'h',h,'file',0);
out.resMCD = zres;
if zraw.objective < obj
    z = zres;
    out.Hsubsets.H1 = zres.Hsubsets.Hopt;
else
    sortmah = sort(mah);
    if h==n
        factor=1;
    else
        factor = sortmah(h)/chi2inv(h/n,k);
    end
    mah = mah/factor;
    weights = mah <= chi2inv(cutoff,k);
    [center_noMCD,cov_noMCD] = weightmecov(X2,weights);
    mah = libra_mahalanobis(X2,center_noMCD,'cov',cov_noMCD);
    z.flag = (mah <= chi2inv(cutoff,k));
    z.center = center_noMCD;
    z.cov = cov_noMCD;
    out.Hsubsets.H1 = is(1:h);
end
covf=z.cov;
centerf=z.center;
[P6,L]=eig(covf);
[L,I]=greatsort(diag(real(L)));
P6=P6(:,I);
out.T=(X2-repmat(centerf,n,1))*P6;
P=P1*P2;
out.P=P(:,1:korig)*P4*P6;
centerfp=center+centerf*rot';
out.M=centerfp;
out.L=L';
out.k=k;
out.h=h;

% creation of Hfreq
out.Hsubsets.Hfreq = zres.Hsubsets.Hfreq(1:h);

outprMCD = out;
%----------------------------------------------------------------------------
function outprAO = projectAO(Tn,k,h,rot,center)

% this function performs the last part of ROBPCA when k is determined.
% input :
%   Tn : the projected data
%   k  : the number of components
%   h  : lower bound for regular observations
%   rot : the rotation matrix
%   center : the classical center of the data

seed=0;
X2=Tn(:,1:k);
[n,p]=size(X2);
outAO=adjustedoutlyingness(X2,'ndir',min(250*p,2500));
AO=outAO.adjout;
cutoffAO=outAO.cutoff;
indexset = find(AO<=cutoffAO)';
%SVD
[P6,T6,L6,r6,Xm6,clmX6]=classSVD(X2(indexset,:));
out.T = (X2-repmat(clmX6,n,1))*P6;
out.P = rot*P6;
out.M = center+clmX6*rot';
out.L = L6;
out.k = k;
out.h = h;
out.AO=AO;
out.cutoff.AO=cutoffAO;
outprAO=out;

%--------------------------------------------------------------------------------
function outDist = CompDist(data,r,out,cutoff,robpcamcd,classic)

% Calculates the distances.
% input: data : the original data
%           r : the rank of the data
%        out is a structure that contains the results of the PCA.
%        classic: an optional structure:
%               classic.P1
%               classic.T1
%               classic.L1
%               classic.clm
%               classic.Xc

if nargin < 6
    options.classic = 0;
else
    options.classic = 1;
end

n = size(data,1);
p = size(data,2);
k = out.k;
skew=out.skew;

% Computing distances
% Robust score distances in robust PCA subspace
if robpcamcd
    if skew==0
        out.sd=sqrt(libra_mahalanobis(out.T,zeros(size(out.T,2),1),'cov',out.L))';
        out.cutoff.sd=sqrt(chi2inv(cutoff,out.k));
    else
        out.sd=out.AO;
        out.cutoff.sd=out.cutoff.AO;
    end
else
    out.sd=zeros(n,1);
    out.cutoff.sd=0;
end
% Orthogonal distances to robust PCA subspace
XRc=data-repmat(out.M,n,1);
Xtilde=out.T*out.P';
Rdiff=XRc-Xtilde;
for i=1:n
    out.od(i,1)=norm(Rdiff(i,:));
end
% Robust cutoff-value for the orthogonal distance
if k~=r
    if skew==0
        [m,s]=unimcd(out.od.^(2/3),out.h);
        out.cutoff.od = sqrt(norminv(cutoff,m,s).^3);
    else
        mcod=mc(out.od);
        if mcod>0
            out.cutoff.od = prctile(out.od,75)+1.5*exp(3*mcod)*iqr(out.od);
        else
            out.cutoff.od = prctile(out.od,75)+1.5*iqr(out.od);
        end
        ttup = sort(-out.od(out.od<out.cutoff.od));
        out.cutoff.od = -ttup(1);
    end
else
    out.cutoff.od=0;
end
if options.classic==1
    % Mahalanobis distance in classical PCA subspace
    Tclas=classic.Xc*classic.P(:,1:out.k);
    out.classic.sd=sqrt(libra_mahalanobis(Tclas,zeros(size(Tclas,2),1),'invcov',1./classic.L(1:out.k)))';
    % Orthogonal distances to classical PCA subspace
    Xtilde=Tclas*classic.P(:,1:out.k)';
    Cdiff=classic.Xc-Xtilde;
    for i=1:n
        out.classic.od(i,1)=norm(Cdiff(i,:));
    end
    out.classic.cutoff.sd=sqrt(chi2inv(cutoff,out.k)); % should be defined after od to have the output in the correct order
    % Classical cutoff-values
    if k~=r
        m=mean(out.classic.od.^(2/3));
        s=sqrt(var(out.classic.od.^(2/3)));
        out.classic.cutoff.od = sqrt(norminv(cutoff,m,s)^3);
    else
        out.classic.cutoff.od=0;
    end
    out.classic.cutoff.sd=sqrt(chi2inv(cutoff,out.k));
    out.classic.flag.od=(out.classic.od<=out.classic.cutoff.od);
    out.classic.flag.sd=(out.classic.sd<=out.classic.cutoff.sd);
    out.classic.flag.all=(out.classic.flag.od)&(out.classic.flag.sd);
    out.classic.class='CPCA';
else
    out.classic=0;
end

out.class='ROBPCA';

if k~=r
    out.flag.od=(out.od<=out.cutoff.od);
    out.flag.sd=(out.sd<=out.cutoff.sd);
    out.flag.all=(out.flag.od)&(out.flag.sd);
else
    out.flag.od=(out.od<=out.cutoff.od);
    out.flag.sd=(out.sd<=out.cutoff.sd);
    out.flag.all=(out.sd<=out.cutoff.sd);
end

outDist = out;

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
%SETPOSTPROC (Re)set POSTPROC field of a datafile
%
%   A = SETPOSTPROC(A,MAPPING)
%   A = SETPOSTPROC(A)
%
% INPUT
%   A        - Datafile
%   POSTPROC - cell containing postprocessing mapping command
%
% OUTPUT
%   A       - Datafile
%
% DESCRIPTION
% Sets the mappings stored in A.POSTPROC. The size of the datafile
% A is set to the output size of PRMAPPING.
% A call without MAPPING clears A.POSTPROC. The size of the datafile
% A is reset to undefined (0).
%
% The mappings in A.POSTPROC may be extended by ADDPOSTPROC.
%
% Mappings in A.POSTPROC are stored only and executed just
% after A is converted from a DATAFILE into a PRDATASET.
%
% SEE ALSO (PRTools Guide)
% DATAFILES, SETPREPROC, ADDPOSTPROC.


text:
function [varargout]=dihedralAngles(varargin)

% function [A]=dihedralAngles(E,V,elementType)
% -----------------------------------------------------------------------
% This function computes the dihedral angles A for the input elements
% defined by E, the element nodal connectivity matrix, and V, the node or
% vertex coordinate matrix. The output array A contains the same number of
% rows as E (a row for each element), and contains m columns where m is the
% number of edges for this type of element (e.g. 12 for a hexahedral element
% or 6 for a tetrahedral element).
%
%
% 2021/03/24 Created KMM
% -----------------------------------------------------------------------

%% PARSE INPUT

switch nargin
    case 2
        E=varargin{1};
        V=varargin{2};
        elementType=[];
    case 3
        E=varargin{1};
        V=varargin{2};
        elementType=varargin{3};
    otherwise
        error('Wrong number of inputs');
end

numNodes=size(E,2);

if isempty(elementType) %have to assume defaults
    switch numNodes
        case 8 %Hexahedral elements
            elementType='hex8';
        case 20 %Hexahedral elements
            elementType='hex20';
            numFaces=6;
        case 4 %Linear tets
            elementType='tet4';
        case 10 %Quadratic tets
            elementType='tet10';
            numFaces=4;
        case 6 %Quadratic triangles
            elementType='tri6';
        case 3
            elementType='tri3';
        case 14
            elementType='rhomdo14';
        otherwise
            elementType='unknown';
    end
end

switch elementType
    case 'tet4'
        numFaces=4;
        c  = [1 3; 3 2; 2 1;...
              1 4; 2 4; 3 4];

        ee = [3 2; 4 3; 1 3;...
              1 2; 4 1; 2 4];
    case 'hex8'
        numFaces=6;
        c=[1 3; 1 4; 1 5; 1 6; ...
           2 3; 2 4; 2 5; 2 6; ...
           3 5; 3 6; 4 5; 4 6; ...
           ];

        ee=[2 1; 4 3; 3 2; 1 4;...
            5 6; 7 8; 6 7; 8 5;...
            2 6; 5 1; 7 3; 4 8;...
            ];
end

%%

[F,~]=element2patch(E,[],elementType);
N=patchNormal(F,V);

FE=reshape(1:1:size(F,1),size(E,1),numFaces);
A=zeros(size(E,1),size(c,1));
for q=1:1:size(c,1)
    a=acos(dot(N(FE(:,c(q,1)),:),N(FE(:,c(q,2)),:),2));
    p1=V(E(:,ee(q,1)),:);
    p2=V(E(:,ee(q,2)),:);
    n=vecnormalize(p2-p1);
    x=cross(N(FE(:,c(q,1)),:),N(FE(:,c(q,2)),:),2);
    logic1=dot(n,x,2)>0;
    logicAngle=(a<pi & logic1);
    a(logicAngle)=pi-a(logicAngle);
    A(:,q)=a;
end

%% Collect output

varargout{1}=A;
if nargout>1
    %Create edge array
    ee=ee';
    ee=ee(:);
    EE=reshape(E(:,ee)',2,numel(ee)*size(E,1)/2)';
    varargout{2}=EE;
end
if nargout>2
    %Create angles for edge array
    AE=reshape(A',1,numel(A))';
    varargout{3}=AE;
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [theta] = normalEqn(X, y)
%NORMALEQN Computes the closed-form solution to linear regression
%   NORMALEQN(X,y) computes the closed-form solution to linear
%   regression using the normal equations.

theta = zeros(size(X, 2), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the code to compute the closed form solution
%               to linear regression and put the result in theta.
%

% ---------------------- Sample Solution ----------------------

theta = pinv((X'*X))*X'*y;

% -------------------------------------------------------------

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

end

author: Ayatans

repo: Machine-Learning-homework


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/ra_mmv_phase_transition_noiseless_s_16.mat';
N = 256;
S = 16;
pta = spx.pursuit.PhaseTransitionAnalysis(N);

% options for CoSaMP MMV solver
solver_options.RankAwareResidual = true;
P = 2;

% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K, S).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K, P, solver_options).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
%isstring_FMAT - Test if parameter is an (admissible) character string.
%
%  USAGE
%
%    test = isstring_FMAT(x,string1,string2,...)
%
%    x              item to test
%    string1...     optional list of admissible strings
%
%  SEE ALSO
%
%    See also isdmatrix, isdvector, isdscalar, isimatrix, isivector, isiscalar.
%

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

function test = isstring_FMAT(x,varargin)

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

test = true;

if ~ischar(x),
	test = false;
	return;
end

if isempty(varargin), return; end

for i = 1:length(varargin),
	if strcmp(x,varargin{i}), return; end
end

test = false;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef prtClassRvmSequential < prtClassRvm
    % prtClassRvmSequential  Relevance vector machine classifier using sequential training
    %
    %   CLASSIFIER = prtClassRvmSequential returns a relevance vector
    %   machine classifier based using sequential training.
    %
    %   CLASSIFIER = prtClassRvmSequential(PROPERTY1, VALUE1, ...)
    %   constructs a prtClassRvmSequential object CLASSIFIER with properties as
    %   specified by PROPERTY/VALUE pairs.
    %
    %    A prtClassRvmSequential object inherits all properties from the abstract class
    %    prtClass. In addition is has the following properties:
    %
    %    kernels                - A cell array of prtKernel objects specifying
    %                             the kernels to use
    %    verbosePlot            - Flag indicating whether or not to plot during
    %                             training
    %    verboseText            - Flag indicating whether or not to output
    %                             verbose updates during training
    %    learningMaxIterations  - The maximum number of iterations
    %
    %    A prtClassRvmSequential also has the following read-only properties:
    %
    %    learningConverged  - Flag indicating if the training converged
    %    beta               - The regression weights, estimated during training
    %    sparseBeta         - The sparse regression weights, estimated during
    %                         training
    %    sparseKernels      - The sparse regression kernels, estimated during
    %                         training
    %
    %   For more information on the algorithm and the above properties, see
    %   the following reference:
    %
    %   Tipping, M. E. and A. C. Faul (2003). Fast marginal likelihood
    %   maximisation for sparse Bayesian models. In C. M. Bishop and B. J.
    %   Frey (Eds.), Proceedings of the Ninth International Workshop on
    %   Artificial Intelligence and Statistics, Key West, FL, Jan 3-6.
    %
    %   prtClassRvmSequential is most useful for datasets with a large
    %   number of observations for which the gram matrix can not be held in
    %   memory. The sequential RVM training algorithm is capable of
    %   operating by generating necessary portions of the gram matrix when
    %   needed. The size of the generated portion of the gram matrix is
    %   determined by the property, largestNumberOfGramColumns. Sequential
    %   RVM training will attempt to generate portions of the gram matrix
    %   that are TraingData.nObservations x largesNumberofGramColums in
    %   size. If the entire gram matrix is this size or smaller it need
    %   only be generated once. Therefore if the entire gram matrix can be
    %   stored in memory, training is much faster. For quickest operation,
    %   largestNumberOfGramColumns should be set as large as possible
    %   without exceeding RAM limitations.
    %
    %    Example:
    %
    %    TestDataSet = prtDataGenUnimodal;      % Create some test and
    %    TrainingDataSet = prtDataGenUnimodal;  % training data classifier
    %    classifier = prtClassRvmSequential('verbosePlot',true); % Create a classifier
    %    classifier = classifier.train(TrainingDataSet);    % Train
    %    classified = run(classifier, TestDataSet);         % Test
    %    % Plot
    %    subplot(2,1,1); classifier.plot;
    %    subplot(2,1,2); [pf,pd] = prtScoreRoc(classified,TestDataSet);
    %    h = plot(pf,pd,'linewidth',3);
    %    title('ROC'); xlabel('Pf'); ylabel('Pd');
    %
    %   See also prtClass, prtClassRvm, prtClassRvnFiguerido,
    %   prtRegressRvmSequential

    properties

        learningPoorlyScaledLikelihoodThreshold = 1e4;
        learningLikelihoodIncreaseThreshold = 1e-6;
        largestNumberOfGramColumns = 5000;
        learningCorrelationRemovalThreshold = 0.99;
        learningFactorRemove = true;   % Remove kernels during train?
        learningRepeatedActionLimit = 25;
    end

    properties(Hidden = true)
        learningResults
    end

    properties (Hidden = true)
        Sigma = [];
    end

    methods
        function Obj = prtClassRvmSequential(varargin)
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end

        function Obj = set.learningPoorlyScaledLikelihoodThreshold(Obj,val)
            assert(prtUtilIsPositiveScalar(val),'prt:prtClassRvmSequential:learningPoorlyScaledLikelihoodThreshold','learningPoorlyScaledLikelihoodThreshold must be a positive scalar');
            Obj.learningPoorlyScaledLikelihoodThreshold = val;
        end
        function Obj = set.learningLikelihoodIncreaseThreshold(Obj,val)
            assert(prtUtilIsPositiveScalar(val),'prt:prtClassRvmSequential:learningLikelihoodIncreaseThreshold','learningLikelihoodIncreaseThreshold must be a positive scalar');
            Obj.learningLikelihoodIncreaseThreshold = val;
        end
        function Obj = set.largestNumberOfGramColumns(Obj,val)
            assert(prtUtilIsPositiveScalarInteger(val),'prt:prtClassRvmSequential:largestNumberOfGramColumns','largestNumberOfGramColumns must be a positive integer scalar');
            Obj.largestNumberOfGramColumns = val;
        end
        function Obj = set.learningCorrelationRemovalThreshold(Obj,val)
            assert(prtUtilIsPositiveScalar(val),'prt:prtClassRvmSequential:learningCorrelationRemovalThreshold','learningCorrelationRemovalThreshold must be a positive scalar');
            Obj.learningCorrelationRemovalThreshold = val;
        end
        function Obj = set.learningFactorRemove(Obj,val)
            assert(prtUtilIsLogicalScalar(val),'prt:prtClassRvmSequential:learningFactorRemove','learningFactorRemove must be a logical scalar');
            Obj.learningFactorRemove = val;
        end
        function Obj = set.learningRepeatedActionLimit(Obj,val)
            assert(prtUtilIsPositiveScalarInteger(val),'prt:prtClassRvmSequential:learningRepeatedActionLimit','learningRepeatedActionLimit must be a positive integer scalar');
            Obj.learningRepeatedActionLimit = val;
        end
    end

    methods (Access=protected, Hidden = true)
        function Obj = trainAction(Obj,DataSet)
            %Rvm = trainAction(Rvm,DataSet) (Private; see prtClass\train)

            warningState = warning;
            %warning off MATLAB:nearlySingularMatrix

            y = Obj.getMinusOneOneTargets(DataSet);

            localKernels = Obj.kernels.train(DataSet);
            nBasis = localKernels.nDimensions;

            if false && nBasis <= Obj.largestNumberOfGramColumns
                Obj = trainActionSequentialInMemory(Obj, DataSet, y);
                return
            end

            if Obj.verboseText
                fprintf('Sequential RVM training with %d possible vectors.\n', nBasis);
            end

            % The sometimes we want y [-1 1] but mostly not
            ym11 = y;
            y(y ==-1) = 0;

            Obj.beta = zeros(nBasis,1);

            relevantIndices = false(nBasis,1); % Nobody!
            alpha = inf(nBasis,1); % Nobody!
            forbidden = zeros(nBasis,1); % Will hold who is forbidding you from joining

            % Find first kernel
            kernelCorrs = zeros(nBasis,1);
            nBlocks = ceil(nBasis./ Obj.largestNumberOfGramColumns);
            for iBlock = 1:nBlocks
                cInds = ((iBlock-1)*Obj.largestNumberOfGramColumns+1):min([iBlock*Obj.largestNumberOfGramColumns nBasis]);

                trainedKernelDownSelected = localKernels.retainKernelDimensions(cInds);
                blockPhi = trainedKernelDownSelected.run_OutputDoubleArray(DataSet);
                blockPhiNormalized = bsxfun(@rdivide,blockPhi,sqrt(sum(blockPhi.*blockPhi))); % We have to normalize here

                kernelCorrs(cInds) = abs(blockPhiNormalized'*ym11);
            end
            [maxVal, maxInd] = max(kernelCorrs); %#ok<ASGLU>

            % Make this ind relevant
            relevantIndices(maxInd) = true;
            selectedInds = maxInd;

            % Add things to forbidden list
            initLogical = false(nBasis,1);
            initLogical(maxInd) = true;
            firstKernel = localKernels.retainKernelDimensions(initLogical);
            firstPhiNormalized = firstKernel.run_OutputDoubleArray(DataSet);
            firstPhiNormalized = firstPhiNormalized - mean(firstPhiNormalized);
            firstPhiNormalized = firstPhiNormalized./sqrt(sum(firstPhiNormalized.^2));
            phiCorrs = zeros(nBasis,1);
            for iBlock = 1:nBlocks
                cInds = ((iBlock-1)*Obj.largestNumberOfGramColumns+1):min([iBlock*Obj.largestNumberOfGramColumns nBasis]);

                trainedKernelDownSelected = localKernels.retainKernelDimensions(cInds);
                if nBlocks > 1 % If there is only one block we can keep the one from before
                    blockPhi = trainedKernelDownSelected.run_OutputDoubleArray(DataSet);
                end
                blockPhiDemeanedNormalized = bsxfun(@minus,blockPhi,mean(blockPhi));
                blockPhiDemeanedNormalized = bsxfun(@rdivide,blockPhiDemeanedNormalized,sqrt(sum(blockPhiDemeanedNormalized.*blockPhiDemeanedNormalized))); % We have to normalize here

                phiCorrs(cInds) = blockPhiDemeanedNormalized'*firstPhiNormalized;
            end
            forbidden(phiCorrs > Obj.learningCorrelationRemovalThreshold) = maxInd;

            % Start the actual Process
            if Obj.verboseText
                fprintf('\t Iteration 0: Intialized with vector %d.\n', maxInd);

                nVectorsStringLength = ceil(log10(length(nBasis)))+1;
            end

            if nBlocks == 1
                % If we have only 1 block we only need to do this once.
                trainedKernelDownSelected = localKernels.retainKernelDimensions(true(nBasis,1));
                PhiM = trainedKernelDownSelected.run_OutputDoubleArray(DataSet);
            end

            % Get the first relevant kernel matrix
            trainedKernelDownSelected = localKernels.retainKernelDimensions(relevantIndices);
            cPhi = trainedKernelDownSelected.run_OutputDoubleArray(DataSet);

            repeatedActionCounter = 0;
            for iteration = 1:Obj.learningMaxIterations

                % Store old log Alpha
                logAlphaOld = log(alpha);

                if iteration == 1
                    % Initial estimates
                    % Estimate Sigma, mu etc.

                    logOut = (ym11*0.9+1)/2;
                    mu = cPhi \ log(logOut./(1-logOut));

                    alpha(relevantIndices) = 1./mu.^2;

                    % Laplacian approx. IRLS
                    A = diag(alpha(relevantIndices));

                    [mu, SigmaInvChol, obsNoiseVar] = prtUtilPenalizedIrls(y,cPhi,mu,A);

                    SigmaChol = inv(SigmaInvChol);
                    Obj.Sigma = SigmaChol*SigmaChol'; %#ok<MINV>

                    yHat = 1 ./ (1+exp(-cPhi*mu));
                end

                % Eval additions and subtractions
                Sm = zeros(nBasis,1);
                Qm = zeros(nBasis,1);
                cError = y-yHat;

                cPhiProduct = bsxfun(@times,cPhi,obsNoiseVar);
                for iBlock = 1:nBlocks
                    cInds = ((iBlock-1)*Obj.largestNumberOfGramColumns+1):min([iBlock*Obj.largestNumberOfGramColumns nBasis]);

                    if nBlocks > 1 % If there is only one block we can keep the on from before
                        trainedKernelDownSelected = localKernels.retainKernelDimensions(cInds);
                        PhiM = trainedKernelDownSelected.run_OutputDoubleArray(DataSet);
                    end

                    Sm(cInds) = (obsNoiseVar'*(PhiM.^2)).' - sum((PhiM.'*cPhiProduct*SigmaChol).^2,2);
                    Qm(cInds) = PhiM.'*cError;
                end

                % Find little sm and qm (these are different for relevant vectors)
                sm = Sm;
                qm = Qm;

                cDenom = (alpha(relevantIndices)-Sm(relevantIndices));
                sm(relevantIndices) = alpha(relevantIndices) .* Sm(relevantIndices) ./ cDenom;
                qm(relevantIndices) = alpha(relevantIndices) .* Qm(relevantIndices) ./ cDenom;

                theta = qm.^2 - sm;
                cantBeRelevent = theta < 0;

                % Addition
                addLogLikelihoodChanges = 0.5*( theta./Sm + log(Sm ./ Qm.^2) ); % Eq (27)
                addLogLikelihoodChanges(cantBeRelevent) = 0; % Can't add things that are disallowed by theta
                addLogLikelihoodChanges(relevantIndices) = 0; % Can't add things already in
                addLogLikelihoodChanges(forbidden > 0) = 0; % Can't add things that are forbidden

                % Removal
                removeLogLikelihoodChanges = -0.5*( qm.^2./(sm + alpha) - log(1 + sm./alpha) ); % Eq (37) (I think this is wrong in the paper. The one in the paper uses Si and Qi, I got this based on si and qi (or S and Q in their code), from corrected from analyzing code from http://www.vectoranomaly.com/downloads/downloads.htm)
                removeLogLikelihoodChanges(~relevantIndices) = 0; % Can't remove things not in
                removeLogLikelihoodChanges(imag(removeLogLikelihoodChanges) > 0) = inf;

                % Modify
                updatedAlpha = sm.^2 ./ theta;
                updatedAlphaDiff = 1./updatedAlpha - 1./alpha;
                modifyLogLikelihoodChanges = 0.5*( updatedAlphaDiff.*(Qm.^2) ./ (updatedAlphaDiff.*Sm + 1) - log(1 + Sm.*updatedAlphaDiff) );

                modifyLogLikelihoodChanges(~relevantIndices) = 0; % Can't modify things not in
                modifyLogLikelihoodChanges(cantBeRelevent) = 0; % Can't modify things that technically shouldn't be in (they would get dropped)

                [addChange, bestAddInd] = max(addLogLikelihoodChanges);
                [remChange, bestRemInd] = max(removeLogLikelihoodChanges);
                [modChange, bestModInd] = max(modifyLogLikelihoodChanges);

                if iteration == 1
                    % On the first iteration we don't allow removal
                    [maxChangeVal, actionInd] = max([addChange, nan, modChange]);
                else
                    if remChange > 0 && Obj.learningFactorRemove
                        % Removing is top priority.
                        % If removing increases the likelihood, we have two
                        % options, actually remove that sample or modify that
                        % sample if that is better
                        [maxChangeVal, actionInd] = max([nan remChange, modifyLogLikelihoodChanges(bestRemInd)]);
                    else
                        % Not going to remove, so we would be allowed to modify
                        [maxChangeVal, actionInd] = max([addChange, remChange, modChange]);
                    end
                end

                if maxChangeVal > Obj.learningPoorlyScaledLikelihoodThreshold
                    warning('prtClassRvm:BadKernelMatrix','Kernel matrix is poorly conditioned. Consider modifying your kernels. Optimization Exiting...' );
                    break
                end

                if maxChangeVal < Obj.learningLikelihoodIncreaseThreshold
                    % There are no good options right now. Therefore we
                    % should exit with the previous iteration stats.
                    Obj.learningConverged = true;
                    Obj.learningResults.exitReason = 'No Good Actions';
                    Obj.learningResults.exitValue = maxChangeVal;
                    if Obj.verboseText
                        fprintf('Convergence criterion met, no necessary actions remaining, maximal change in log-likelihood %g\n\n',maxChangeVal);
                    end

                    break;
                end

                switch actionInd
                    case 1
                        cBestInd = bestAddInd;
                    case 2
                        cBestInd = bestRemInd;
                    case 3
                        cBestInd = bestModInd;
                end
                if iteration > 1 && lastAction == actionInd && cBestInd == lastInd
                    repeatedActionCounter = repeatedActionCounter + 1;
                else
                    repeatedActionCounter = 0;
                end

                if repeatedActionCounter >= Obj.learningRepeatedActionLimit
                    if Obj.verboseText
                        fprintf('Exiting... repeating action limit has been reached.\n\n');
                    end
                    return
                end

                if Obj.verboseText
                    actionStrings = {sprintf('Addition: Vector %s has been added.  ', sprintf(sprintf('%%%dd',nVectorsStringLength),bestAddInd));
                                     sprintf('Removal:  Vector %s has been removed.', sprintf(sprintf('%%%dd',nVectorsStringLength), bestRemInd));
                                     sprintf('Update:   Vector %s has been updated.', sprintf(sprintf('%%%dd',nVectorsStringLength), bestModInd));};
                    fprintf('\t Iteration %d: %s Change in log-likelihood %g.\n',iteration, actionStrings{actionInd}, maxChangeVal);
                end

                lastAction = actionInd;
                switch actionInd
                    case 1 % Add

                        relevantIndices(bestAddInd) = true;
                        selectedInds = cat(1,selectedInds,bestAddInd);

                        alpha(bestAddInd) = updatedAlpha(bestAddInd);
                        % Modify Mu
                        % (Penalized IRLS will fix it soon but we need good initialization)

                        trainedKernelDownSelected = localKernels.retainKernelDimensions(bestAddInd);
                        newPhi = trainedKernelDownSelected.run_OutputDoubleArray(DataSet);

                        cFactor = (Obj.Sigma*(cPhi)'*(newPhi.*obsNoiseVar));
                        Sigmaii = 1./(updatedAlpha(bestAddInd) + Sm(bestAddInd));
                        newMu = Sigmaii*Qm(bestAddInd);

                        updatedOldMu = mu - newMu*cFactor;
                        sortedSelected = sort(selectedInds);
                        newMuLocation = find(sortedSelected==bestAddInd);

                        mu = zeros(length(mu)+1,1);
                        mu(setdiff(1:length(mu),newMuLocation)) = updatedOldMu;
                        mu(newMuLocation) = newMu;

                        % Add things to forbidden list
                        newPhiDemeanedNormalized = newPhi - mean(newPhi);
                        newPhiDemeanedNormalized = newPhiDemeanedNormalized./sqrt(sum(newPhiDemeanedNormalized.^2));
                        phiCorrs = zeros(nBasis,1);
                        for iBlock = 1:nBlocks
                            cInds = ((iBlock-1)*Obj.largestNumberOfGramColumns+1):min([iBlock*Obj.largestNumberOfGramColumns nBasis]);

                            if nBlocks > 1
                                trainedKernelDownSelected = localKernels.retainKernelDimensions(cInds);
                                blockPhiDemeanedNormalized = trainedKernelDownSelected.run_OutputDoubleArray(DataSet);

                                blockPhiDemeanedNormalized = bsxfun(@minus,blockPhiDemeanedNormalized,mean(blockPhiDemeanedNormalized));
                                blockPhiDemeanedNormalized = bsxfun(@rdivide,blockPhiDemeanedNormalized,sqrt(sum(blockPhiDemeanedNormalized.*blockPhiDemeanedNormalized))); % We have to normalize here

                            %else we have this from before
                            end

                            phiCorrs(cInds) = blockPhiDemeanedNormalized'*newPhiDemeanedNormalized;
                        end
                        forbidden(phiCorrs > Obj.learningCorrelationRemovalThreshold) = bestAddInd;

                        lastInd = bestAddInd;
                    case 2 % Remove

                        removingInd = sort(selectedInds)==bestRemInd;

                        mu = mu + mu(removingInd) .* Obj.Sigma(:,removingInd) ./ Obj.Sigma(removingInd,removingInd);
                        mu(removingInd) = [];

                        relevantIndices(bestRemInd) = false;
                        selectedInds(selectedInds==bestRemInd) = [];
                        alpha(bestRemInd) = inf;

                        % Anything this guy said is forbidden is now
                        % allowed.
                        forbidden(forbidden == bestRemInd) = 0;

                        lastInd = bestRemInd;

                    case 3 % Modify
                        modifyInd = sort(selectedInds)==bestModInd;

                        alphaChangeInv = 1/(updatedAlpha(bestModInd) - alpha(bestModInd));
                        kappa = 1/(Obj.Sigma(modifyInd,modifyInd) + alphaChangeInv);
                        mu = mu - mu(modifyInd)*kappa*Obj.Sigma(:,modifyInd);

                        alpha(bestModInd) = updatedAlpha(bestModInd);

                        lastInd = bestModInd;
                end

                % At this point relevantIndices and alpha have changes.
                % Now we re-estimate Sigma, mu, and sigma2
                if nBlocks > 1
                    trainedKernelDownSelected = localKernels.retainKernelDimensions(relevantIndices);
                    cPhi = trainedKernelDownSelected.run_OutputDoubleArray(DataSet);
                else
                    cPhi = PhiM(:,relevantIndices);
                end

                % Laplacian approx. IRLS
                A = diag(alpha(relevantIndices));

                if isempty(cPhi)
                    yHat = 0.5*ones(size(y));
                else
                    [mu, SigmaInvChol, obsNoiseVar] = prtUtilPenalizedIrls(y,cPhi,mu,A);

                    SigmaChol = inv(SigmaInvChol);
                    Obj.Sigma = SigmaChol*SigmaChol'; %#ok<MINV>

                    yHat = 1 ./ (1+exp(-cPhi*mu)); % We use a logistic here.
                end

                % Store beta
                Obj.beta = zeros(nBasis,1);
                Obj.beta(relevantIndices) = mu;

                if ~mod(iteration,Obj.verbosePlot)
                    if DataSet.nFeatures == 2
                        Obj.verboseIterationPlot(DataSet,relevantIndices);
                    elseif iteration == 1
                        warning('prt:prtClassRvmSequential','Learning iteration plot can only be produced for training Datasets with 2 features');
                    end
                end

                % Check tolerance
                changeVal = abs(log(alpha)-logAlphaOld);
                changeVal(isnan(changeVal)) = 0; % inf-inf = nan
                if all(changeVal < Obj.learningConvergedTolerance) && iteration > 1
                    Obj.learningConverged = true;
                    Obj.learningResults.exitReason = 'Alpha Not Changing';
                    Obj.learningResults.exitValue = TOL;
                    if Obj.verboseText
                        fprintf('Exiting...Precisions no longer changing appreciably.\n\n');
                    end
                    break;
                end
            end

            if Obj.verboseText && iteration == Obj.learningMaxIterations
                fprintf('Exiting...Convergence not reached before the maximum allowed iterations was reached.\n\n');
            end

            % Make sparse represenation
            Obj.sparseBeta = Obj.beta(relevantIndices,1);
            Obj.sparseKernels = localKernels.retainKernelDimensions(relevantIndices);

            % Very bad training
            if isempty(Obj.sparseBeta)
                warning('prt:prtClassRvm:NoRelevantFeatures','No relevant features were found during training.');
            end

            % Reset warning
            warning(warningState);

        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
if(exist('D:\Datasets\fera/au_training', 'file'))
    FERA2011_dir = 'D:\Datasets\fera/au_training/';
elseif(exist('/multicomp/datasets/fera2011/', 'file'))
    FERA2011_dir = '/multicomp/datasets/fera2011/au_training/';
elseif(exist('E:\Datasets\fera/au_training', 'file'))
    FERA2011_dir = 'E:\Datasets\fera/au_training/';
elseif(exist('/media/tadas/2EBEA130BEA0F20F/datasets/fera/au_training/', 'file'))
    FERA2011_dir = '/media/tadas/2EBEA130BEA0F20F/datasets/fera/au_training/';
else
    fprintf('FERA2011 location not found (or not defined)\n');
end

all_recs = {'train_001', 'train_002', 'train_003', 'train_004', 'train_005',...
            'train_006', 'train_007', 'train_008', 'train_009', 'train_010',...
            'train_011', 'train_012', 'train_013', 'train_014', 'train_015',...
            'train_016', 'train_017', 'train_018', 'train_019', 'train_020',...
            'train_021', 'train_022', 'train_023', 'train_024', 'train_025',...
            'train_026', 'train_027', 'train_028', 'train_029', 'train_030',...
            'train_031', 'train_032', 'train_033', 'train_034', 'train_035',...
            'train_036', 'train_037', 'train_038', 'train_039', 'train_040',...
            'train_041', 'train_042', 'train_043', 'train_044', 'train_045',...
            'train_046', 'train_047', 'train_048', 'train_049', 'train_050',...
            'train_051', 'train_052', 'train_053', 'train_054', 'train_055',...
            'train_056', 'train_057', 'train_058', 'train_059', 'train_060',...
            'train_061', 'train_062', 'train_063', 'train_064', 'train_065',...
            'train_066', 'train_067', 'train_068', 'train_069', 'train_070',...
            'train_071', 'train_072', 'train_073', 'train_074', 'train_075',...
            'train_076', 'train_077', 'train_078', 'train_079', 'train_080',...
            'train_081', 'train_082', 'train_083', 'train_084', 'train_085',...
            'train_086', 'train_087'};

% Making them person independent
train_recs = all_recs(1:56);
devel_recs = setdiff(all_recs, train_recs);

all_aus = [1, 2, 4, 6, 7, 10, 12, 15, 17, 18, 25, 26];

author: TadasBaltrusaitis

repo: OpenFace


text:
function stroud_test085 ( )

%*****************************************************************************80
%
%% TEST085 tests CIRCLE_ANNULUS_AREA_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  ntest = 3;

  radius1_test = [ 0.0, 1.0, 1.0 ];
  radius2_test = [ 1.0, 2.0, 3.0 ];
  xc_test = [ 0.0, 1.0, 3.0 ];
  yc_test = [ 0.0, 0.0, 4.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST085\n' );
  fprintf ( 1, '  CIRCLE_ANNULUS_AREA_2D computes the area of a \n' );
  fprintf ( 1, '    circular annulus.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     XC        YC          Radius1    Radius2    Area\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest

    xc = xc_test(i);
    yc = yc_test(i);

    radius1 = radius1_test(i);
    radius2 = radius2_test(i);

    area = circle_annulus_area_2d ( radius1, radius2 );

    fprintf ( 1, '  %9f  %9f  %9f  %9f  %9f\n', xc, yc, radius1, radius2, area );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Az, El, D] = topocent(X,dx)
%TOPOCENT  Transformation of vector dx into topocentric coordinate
%          system with origin at X.
%          Both parameters are 3 by 1 vectors.
%          Output: D    vector length in units like the input
%                  Az   azimuth from north positive clockwise, degrees
%                  El   elevation angle, degrees

%Kai Borre 11-24-96
%Copyright (c) by Kai Borre
%$Revision: 1.0 $  $Date: 1997/09/26  $

dtr = pi/180;
[phi,lambda,~] = togeod(6378137,298.257223563,X(1),X(2),X(3));
cl = cos(lambda*dtr); sl = sin(lambda*dtr);
cb = cos(phi*dtr); sb = sin(phi*dtr);
F = [-sl -sb*cl cb*cl;
      cl -sb*sl cb*sl;
       0    cb   sb];
local_vector = F'*dx;
E = local_vector(1);
N = local_vector(2);
U = local_vector(3);
hor_dis = sqrt(E^2+N^2);
if hor_dis < 1.e-20
   Az = 0;
   El = 90;
else
   Az = atan2(E,N)/dtr;
   El = atan2(U,hor_dis)/dtr;
end
if Az < 0
   Az = Az+360;
end
D = sqrt(dx(1)^2+dx(2)^2+dx(3)^2);
%%%%%%%%% end topocent.m %%%%%%%%%

author: yandld


text:
function data = montage_image_Worsley(image_name, varargin)
% ::
%
%    data = montage_image_Worsley(3D or 4D image name)
%
% Make a compact montage of some images
% Designed by Keith Worsley for pca_image.m
% Adapted by Tor Wager, Feb 2008
%
% Limits and Colormap are designed for component loadings between [-1 1]
%
% :Examples:
% ::
%
%    create_figure('Montage');
%    montage_image_Worsley('test_run1_pca.img');
%
%    data = montage_image_Worsley(imgname, 'pcacov')  % changes scaling and colormap
%    data = montage_image_Worsley(imgname, 'pcacov', [1 3 5])  % show only
%                                                           volumes 1, 3, 5 in image(s)

% ..
%    COPYRIGHT:   Copyright 2002 K.J. Worsley,
%              Department of Mathematics and Statistics,
%              McConnell Brain Imaging Center,
%              Montreal Neurological Institute,
%              McGill University, Montreal, Quebec, Canada.
%              worsley@math.mcgill.ca
%
%              Permission to use, copy, modify, and distribute this
%              software and its documentation for any purpose and without
%              fee is hereby granted, provided that this copyright
%              notice appears in all copies. The author and McGill University
%              make no representations about the suitability of this
%              software for any purpose.  It is provided "as is" without
%              express or implied warranty.
% ..
%
% ..
%    Tor Wager
% ..

scaletype = 'pcacor';  % Scale image and colors for PCA component viewing
whichvols = [];

for i = 1:length(varargin)
        if ischar(varargin{i})
            switch varargin{i}
                % reserved keywords
                case 'noscale', scaletype = 'none';

                case 'pcacov', scaletype = 'pcacov';

                case 'whichvols', whichvols = varargin{i + 1};

                otherwise, warning(['Unknown input string option:' varargin{i}]);
            end
        end
end

% Get data
% -------------------------------------------------------------------------
image_name = expand_4d_filenames(image_name); % for SPM2 compatibility
V = spm_vol(image_name);

if isempty(whichvols)
    wh_vols = 1:length(V);

else
    % select only certain vols
    V = V(wh_vols);
end

p = length(V);
go_ok = 1;
if p > 10
    go_ok = input('More than 10 images! Are you sure you want to continue?');
end
if ~go_ok, return, end

data = spm_read_vols(V);

% Setup stuff
% -------------------------------------------------------------------------

numslices = V(1).dim(3);
% p = num comps
numys=V(1).dim(2);
numxs=V(1).dim(1);
   xr=1:numxs;
   yr=1:numys;
   %N=prod(V(1).dim(1:3));

nrow=round(sqrt(numslices/3.25/p));
nrow=max(nrow,1);

ncol=ceil(numslices/nrow);

numxr=length(xr);
numyr=length(yr);
bigmat=zeros(p*nrow*numyr,ncol*numxr);

% Create big mat
% -------------------------------------------------------------------------
r=0;
for k=1:p
   c=0;
   for i=1:numslices
      bigmat((1:numyr)+r*numyr,(1:numxr)+c*numxr)=flipud(data(xr,yr,i,k)');
      c=c+1;
      if c>=ncol && i<numslices
         r=r+1;
         c=0;
      end
   end
   r=r+1;
end
xtick=((1:ncol*numxr)-1)/(ncol*numxr-1)*ncol-0.5;
ytick=((1:p*nrow*numyr)-1)/(p*nrow*numyr-1)*p+0.5;

switch scaletype
    case 'pcacor'
        zlim=[-1 1]; % [min(bigmat(:)) max(bigmat(:))];
    case 'pcacov'
        sdbm = nanstd(data(:));
        if(sdbm == 0)
            warning('Std dev of 0!');
            zlim = [-3 3];
        else
            zlim = [-3*sdbm 3*sdbm];
        end

    otherwise
        zlim = [min(bigmat(:)) max(bigmat(:))];
end

%
% Make image
% -------------------------------------------------------------------------
imagesc(xtick,ytick,bigmat,zlim);
xlabel('Slice (0 based)'); ylabel('Component');
switch scaletype
    case {'pcacor', 'pcacov'}
        title('Spatial components');
    otherwise
end
%colormap(spectral); colorbar;

set(gca,'YDir','Reverse')

% tor colormap stuff
% -------------------------------------------------------------------------
switch scaletype
    case {'pcacor'}
    cmap = colormap_tor([0 0 0], [1 1 0], [0 0 1], [0 1 1], [.5 .5 .5], [1 .4 .4], [1 0 0]);
    case {'pcacov'}
       cmap = colormap_tor([0 0 .5], [1 1 0], [0 0 1], [0 1 1], [0 0 0], [1 .4 .4], [1 0 0]);

    otherwise
     cmap = colormap_tor([0 0 0], [1 1 0], [0 0 1], [0 1 1], [.5 .5 .5], [1 .4 .4], [1 0 0]);
end
colormap(cmap)
colorbar('vert')

axis tight

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function HamilExp = DoublePendulumEnergy(in1,in2)
%DOUBLEPENDULUMENERGY
%    HAMILEXP = DOUBLEPENDULUMENERGY(IN1,IN2)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    21-Sep-2019 22:29:59

dz1 = in2(:,1);
dz2 = in2(:,2);
z1 = in1(:,1);
z2 = in1(:,2);
HamilExp = cos(z1).*(2.09e2./2.0e2)+cos(z2).*3.268e-1+dz1.^2.*1.376e-2+dz2.^2.*3.272e-3+dz1.*dz2.*cos(z1-z2).*8.885000000000001e-3;

author: dynamicslab

repo: SINDy-PI


text:
function r = sinvchi2rand(nu, s2, M, N)
% SINVCHI2RAND  Random matrices from scaled inverse-chi distribution
%
%  R = SINVCHI2RAND(NU, S2)
%  R = SINVCHI2RAND(NU, S2, M, N)
%
%  Returns a randon number/matrix R from scaled inverse-chi square
%  distribution. Nu is the degrees of freedom and S2 is the scale
%  squared. Parametrisation is according to Gelman et. al. (2004).

% Copyright (c) 1998-2004 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 2
  error('Too few arguments');
end
if nargin==2
    [M,N]=size(s2);
else
    if numel(s2)>1 || numel(nu)>1
        error('Arguments M and N can only be used if nu and s2 are scalars');
    end
    if nargin < 4
        N=1;
    end
end
r=nu.*s2./chi2rnd(nu,M,N);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function out = saveobj(obj)
%SAVEOBJ (overloaded)
warning('YALMIP objects cannot be saved in binary format. You will run into troubles if you try to load this file later. You should clear all YALMIP objects first, or avoid having any YALMIP objects in the list of variables which you save.');
out = [];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef Quadrature_Quadrilateral < Quadrature

    methods (Access = public)
        function computeQuadrature(obj,order)
            computeQuadrature@Quadrature(obj,order);
            switch order
                case 'CONSTANT'
                    obj.ngaus = 1;
                    obj.posgp(:,1) = [0,0];
                    obj.weigp = 4;

                case 'LINEAR'
                    obj.ngaus = 4;
                    % Compute WEIGP and POSGP
                    a =  0.577350269189626;
                    obj.posgp(:,1) = [-a,-a];
                    obj.posgp(:,2) = [+a,-a];
                    obj.posgp(:,3) = [-a,+a];
                    obj.posgp(:,4) = [+a,+a];

                    obj.weigp =  [1,1,1,1];%1*ones(1,obj.ngaus);

                case 'QUADRATIC' %SERENDIPITY, QUADRILATERAL QUADRATIC NOT IMPLEMENTED
                    % Copied from down below
                    obj.ngaus = 9;
                    % Compute WEIGP and POSGP
                    a =  0.77459667;
                    obj.posgp(:,1) = [ 0,+a];
                    obj.posgp(:,2) = [ 0, 0];
                    obj.posgp(:,3) = [+a,+a];
                    obj.posgp(:,4) = [-a,-a];
                    obj.posgp(:,5) = [-a, 0];
                    obj.posgp(:,6) = [+a, 0];
                    obj.posgp(:,7) = [+a,-a];
                    obj.posgp(:,8) = [-a,+a];
                    obj.posgp(:,9) = [ 0,-a];

                    obj.weigp =(4/9)*ones(1,obj.ngaus);

                case 'QUADRATICMSS' %SERENDIPITY, QUADRILATERAL QUADRATIC NOT IMPLEMENTED
                    obj.ngaus = 9;
                    % Compute WEIGP and POSGP
                    a =  0.77459667;
                    obj.posgp(:,1) = [ 0,+a];
                    obj.posgp(:,2) = [ 0, 0];
                    obj.posgp(:,3) = [+a,+a];
                    obj.posgp(:,4) = [-a,-a];
                    obj.posgp(:,5) = [-a, 0];
                    obj.posgp(:,6) = [+a, 0];
                    obj.posgp(:,7) = [+a,-a];
                    obj.posgp(:,8) = [-a,+a];
                    obj.posgp(:,9) = [ 0,-a];

                    obj.weigp =(4/9)*ones(1,obj.ngaus);

                case 'QUADRATICMASS'
                    posgl(1) =-0.774596669241483;
                    posgl(2) = 0.0;
                    posgl(3) = 0.774596669241483;
                    weigl(1) = 0.555555555555556;
                    weigl(2) = 0.888888888888889;
                    weigl(3) = 0.555555555555556;
                    obj.ngaus = 9;
                    igaus = 0;
                    nlocs = 3;
                    for ilocs = 1:nlocs
                        for jlocs = 1:nlocs
                            igaus = igaus+1;
                            obj.weigp(  igaus) = weigl(ilocs)*weigl(jlocs);
                            obj.posgp(1,igaus) = posgl(ilocs);
                            obj.posgp(2,igaus) = posgl(jlocs);
                        end
                    end

                otherwise
                    error('Invalid interpolation order for element Quadrilateral.');
            end
        end
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [oz] = l2oz(l)
% Convert volume from liters to US liquid ounces.
% Chad Greene 2012
oz = l*33.814022701;

author: Sable

repo: mcbench-benchmarks


text:
function calpak_test622 ( )

%*****************************************************************************80
%
%% CALPAK_TEST622 tests YEAR_TO_DOMINICAL_JULIAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST622\n' );
  fprintf ( 1, '  For the Julian calendar,\n' );
  fprintf ( 1, '  YEAR_TO_DOMINICAL_JULIAN determines the dominical number.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Year  Dominical Number\n' );
  fprintf ( 1, '\n' );

  for y = 1577 : 1587

    s = y_to_s ( y );
    [ n1, n2 ] = year_to_dominical_julian ( y );
    d1 = i4_to_a ( n1 );

    if ( n1 == n2 )
      fprintf ( 1, '  %10s  %d  %s\n', s, n1, d1 );
    else
      d2 = i4_to_a ( n2 );
      fprintf ( 1, '  %10s  %d  %s  %d  %s\n', s, n1, d1, n2, d2 );
    end

  end

  return
end


text:
function cmp = yjf_compare ( y1, j1, f1, y2, j2, f2 )

%*****************************************************************************80
%
%% YJF_COMPARE compares two YJF dates.
%
%  Discussion:
%
%    The routine is "generic" and does not assume a particular calendar.
%    However, it does assume that the calendar dates are "normalized".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, J1, real F1,
%    the first YJF date.
%
%    Input, integer Y2, J2, real F2,
%    the second YJF date.
%
%    Output, character CMP:
%    '<' if date 1 precedes date 2;
%    '=' if date 1 equals date 2;
%    '>' if date 1 follows date 2;
%
  if ( y1 < y2 )
    cmp = '<';
  elseif ( y1 > y2 )
    cmp = '>';
  else

    if ( j1 < j2 )
      cmp = '<';
    elseif ( j1 > j2 )
      cmp = '>';
    else

      if ( f1 < f2 )
        cmp = '<';
      elseif ( f1 > f2 )
        cmp = '>';
      else
        cmp = '=';
      end

    end

  end

  return
end


text:
function [oz] = mi32oz(mi3)
% Convert volume from cubic miles to US liquid ounces.
% Chad Greene 2012
oz = mi3*140942994870000;