text:
classdef VS_Stage < VS_AbstractStage
    %VS_Stage Basic stage with one tank and one engine (and no boosters)
    %   Detailed explanation goes here

    properties
        %direct inputs
        isp double = 300; %seconds
        dryMassFrac(1,1) double = 0.2;
        desiredT2W(1,1) double = 1.0;
        bodyInfo KSPTOT_BodyInfo = KSPTOT_BodyInfo.empty(1,0);

        %computed
        %payloadStage(1,:) VS_Stage = VS_Stage.empty(1,0);

        %optimized quantities
        dryMass(1,1) double = 1; %kg
        deltaV(1,1) double = 0; %m/s

        %flags
        isPayload(1,1) logical = false;
    end

    methods
        function obj = VS_Stage(phase)
            obj.phase = phase;
        end

        function deltaV = getStageDeltaV(obj)
            deltaV = obj.deltaV;
        end

        function dryMass = getTotalStageDryMass(obj)
            dryMass = obj.dryMass;
        end

        function expPropMass = computeExpendedPropMass(obj)
            expPropMass = obj.computeInitMass() - obj.computeFinalMass();
        end

        function tgtDryMass = computeTgtDryMass(obj)
            tgtDryMass = (obj.dryMassFrac/(1 - obj.dryMassFrac)) * obj.computeExpendedPropMass();
        end

        function stageTotalMass = computeStageTotalMass(obj)
            stageTotalMass = obj.dryMass + obj.computeExpendedPropMass();
        end

        function str = getStageOutputStr(obj)
            str = cell(0,1);
            str{end+1} = sprintf('\tStage: %s', obj.name);
            str{end+1} = sprintf('\t\tTotal Stage Dry Mass: %0.3f mT', obj.dryMass/1000);
            str{end+1} = sprintf('\t\tTotal Stage Prop Mass: %0.3f mT', obj.computeExpendedPropMass()/1000);
            str{end+1} = sprintf('\t\tTotal Stage Mass: %0.3f mT', obj.computeStageTotalMass()/1000);
            str{end+1} = sprintf('\t\tTotal Stage Delta-V: %0.3f km/s', obj.deltaV/1000);
            str{end+1} = sprintf('\t\tStage Req''d Thrust: %0.3f kN', obj.computeReqdThrust()/1000);
            str{end+1} = sprintf('\t\tStage Burn Time: %0.3f sec', obj.computeEngineBurnTime());
            str{end+1} = sprintf('\t\t----');
            str{end+1} = sprintf('\t\tStage Isp: %0.3f sec', obj.isp);
            str{end+1} = sprintf('\t\tStage Dry Mass Fraction: %0.3f', obj.dryMassFrac);
            str{end+1} = sprintf('\t\tStage Celestial Body: %s', obj.bodyInfo.name);
        end

        %%%%%%%% For Optimizer Use %%%%%%%%%
        function initStageForOpt(obj, phaseDvPerStage)
            obj.deltaV = phaseDvPerStage;

            if(not(obj.isPayloadOnlyStage()))
                diff = Inf;
                iterCnt = 1;
                while(diff >= 0.01 && iterCnt <=200)
                    obj.dryMass = obj.computeTgtDryMass();

                    diff = abs(obj.dryMass - obj.computeTgtDryMass());
                    iterCnt = iterCnt + 1;
                end
            else
                if(obj.dryMassFrac == 1)
                    obj.deltaV = 0;
                end
            end
        end

        function x = getStageOptimizerVars(obj)
            x = [];

            if(not(obj.isPayloadOnlyStage()))
                x = [obj.deltaV, obj.dryMass];
            else
                x = [obj.deltaV];
            end
        end

        function updateStageWithVars(obj, x)
            if(not(obj.isPayloadOnlyStage()))
                obj.deltaV = x(1);
                obj.dryMass = x(2);
            else
                obj.deltaV = x(1);
            end
        end

        function [lb,ub] = getStageOptVarBounds(obj)
            lb = [];
            ub = [];

            if(not(obj.isPayloadOnlyStage()))
                lb = [0, 0]; %dV and dry mass have lower bounds of 0
                ub = [Inf, Inf]; %no upper bounds on these variables
            else
                lb = 0;
                ub = Inf;
            end
        end

        function [c, ceq] = getStageNlConValues(obj)
            c = [];
            ceq = [];

%             if(not(obj.isPayloadOnlyStage()))
                ceq(1) = obj.getTotalStageDryMass() - obj.computeTgtDryMass();
%             end
        end
    end

    methods(Access=private)
        function finalMass = computeFinalMass(obj)
            finalMass = obj.computePayloadMass() + obj.dryMass;
        end

        function initMass = computeInitMass(obj)
            finalMass = obj.computeFinalMass();
            g0 = getG0();

            initMass = finalMass * exp(obj.deltaV/(g0*obj.isp));
        end

        function payloadMass = computePayloadMass(obj)
            if(isempty(obj.payloadStage))
                payloadMass = 0;
            else
                payloadMass = obj.payloadStage.computeInitMass();
            end
        end

        function reqdThrust = computeReqdThrust(obj)
            reqdThrust = obj.desiredT2W * (obj.computeStageTotalMass() + obj.computePayloadMass()) * obj.computeGforBody();
        end

        function engineBurnTime = computeEngineBurnTime(obj)
            engineBurnTime = obj.computeExpendedPropMass() / (obj.computeReqdThrust() / (getG0() * obj.isp));
        end

        function tf = isPayloadOnlyStage(obj)
            tf = obj.isPayload;
        end
    end
end


text:
% SHUFFLE - shuffle a given dimension in an array
%
% Usage: >> Y = shuffle(X)
%        >> [Y = shuffle(X, DIM)
%
% Inputs:
%   X   - input array
%   DIM - dimension index (default is first non-singleton dimension)
%
% Outputs:
%    Y - shuffled array
%    I - forward indices (Y = X(I) if 1D)
%    J - reverse indices (X(J) = Y if 1D)
%
% Author: Arnaud Delorme, SCCN/INC/UCSD USA, Dec 2000

% Copyright (C) Arnaud Delorme, SCCN/INC/UCSD USA, Dec 2000
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [x, i, j]=shuffle( y, dim)

if nargin < 1
	help shuffle;
	return;
end
if nargin < 2
	if size(y,1) ~= 1
		dim = 1;
	else
		if size(y,2) ~= 1
			dim = 2;
		else
			dim = 3;
		end
	end
end

r =size(y, dim);
a = rand(1,r);
[tmp i] = sort(a);
switch dim
	case 1
		x = y(i,:,:,:,:);
	case 2
		x = y(:,i,:,:,:);
	case 3
		x = y(:,:,i,:,:);
	case 4
		x = y(:,:,:,i,:);
	case 5
		x = y(:,:,:,:,i);
end;
[tmp j] = sort(i); % unshuffle

return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [ spiketimes_jitt ] = JitterSpiketimes(spiketimes,jitterwin,varargin)
%[ spiketimes ] = JitterSpiketimes(spiketimes,jitterwin) jitters spiketimes
%of all cells within a designated window.
%
%INPUT
%   spiketimes  {Ncells} cell array of [Nspikes] vectors of spiketimes for
%               each cell (seconds)
%   jitterwin   time window within which to jitter spike times (seconds)
%
%   (options)
%       'jittertype'    'window' or 'centered'.
%                       -'centered' jitters each spike around its time.
%                       -'window' jitters within its window (more rigorous, default)
%                           see Jonathan Platkewitz paper
%
%OUTPUT
%   spiketimes  jittered spiketimes
%
%
%DLevenstein 2016
%%
%parse args
p = inputParser;
addParameter(p,'jittertype','window')

parse(p,varargin{:})
jittertype = p.Results.jittertype;

%%
numcells = length(spiketimes);
if isa(spiketimes,'tsdArray')
    for c = 1:numcells
        spiketimestemp{c} = Range(spiketimes{c},'s');
    end
    spiketimes = spiketimestemp;
    clear spiketimestemp
end

%Do this with cellfun....
switch jittertype
    case 'centered'
        for c = 1:numcells
            spiketimes_jitt{c} = spiketimes{c}+2*jitterwin*rand(size(spiketimes{c}))-jitterwin;
        end
    case 'window'
        %Question: add global offset, so not always starting windows at same
        %points in recording?
        spiketimes_jitt = cellfun(@(X) ...
            jitterwin.*(floor(X/jitterwin)+rand(size(X))),...
            spiketimes,'UniformOutput',false);
end

%Re-sort
spiketimes_jitt = cellfun(@(X) sort(X),spiketimes_jitt,'UniformOutput',false);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function s=pval_to_star(p)
if p<0.001
    s='***';
elseif p<0.01
    s='**';
elseif p<0.05
    s='*';
else
    s='';
end
end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
function makeplot(out,varargin)

%MAKEPLOT makes plots for the main functions. These figures can also be obtained
% by setting 'plots = 1' in those functions.
%
% Required input:
%  out = a structure containing the output of one of the following classes:
%        MCDCOV, LS, LTS, MLR, MCDREG, CPCA,CPCR, CSIMPLS, RAPCA, ROBPCA,
%        RPCR, RSIMPLS, CDA, RDA
%
% Optional input:
%  nameplot  :  0         : menu of plots     (default = 0)
%              'all'      : all possible plots
%              'scree'    : scree plot
%              'pcadiag'  : Score outlier map
%              'regdiag'  : Regression outlier map
%              '3ddiag'   : 3D Regression outlier map
%              'robdist'  : A plot of the robust distances
%              'qqmcd'    : A qq-plot of the robust distances versus the quantiles of the chi-squared distribution
%              'dd'       : A DD-plot: Robust distances versus Mahalanobis distances
%              'ellipse'  : A tolerance ellipse
%              'resfit'   : Standardized LTS Residuals versus fitted values
%              'resindex' : Standardized LTS Residuals versus index
%              'qqlts'    : Normal QQ-plot of the LTS residuals
%              'scatter'  : Scatterplot with LTS line
%              'da'       : Tolerance ellipses of a discrimant analysis
%              'simca'    : Scatterplot with boundaries defined by the number of PC's from a simca method
%  labod   : number of points to be identified in score plots
%            with largest orthogonal distance.
%  labsd   : number of points to be identified in score and regression plots
%            with largest score distance.
%  labresd : number of points to be identified in regression plots
%            with largest residual distance.
%  labmcd  : number of points to be identified in MCD plots
%            with the largest robust distance.
%  lablts  : number of points to be identified in LTS plots
%            with the largest absolute standardized residual.
%  classic : In case of a robust analysis, classic can be set to 0 to avoid classical plots (default).
%            If set to one, the inputargument 'out' must contain a field called 'classic' which is a structure.
%
% I/O: makeplot(out,'nameplot',0,'labsd',3,'labod',3,'labresd',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.
%
% Example: outrpcr=rpcr(x,y,'plots',0,'classic',1);
%          makeplot(outrpcr,'labsd',5, 'classic',1)
%          makeplot(outrpcr,'nameplot','scree')
%
%          outmcd=mcdcov(x);
%          makeplot(outmcd,'labmcd',6)
%          makeplot(outmcd,'labmcd',4,'nameplot','dd')
%
%          outls=libra_ols(x,y);
%          makeplot(outls)
%
% This function is part of the Matlab Library for Robust Analysis (LIBRA),
% available at:
%              www.wis.kuleuven.ac.be/stat/robust.html
%
% Written by Sabine Verboven  09/04/2003
% Last Updated : 29/01/2008
%
% Uses functions:  screeplot, scorediagplot, regresdiagplot,
%                  regresdiagplot3d, plotnumbers, plotnumbers3d,
%                  distplot, chiqqplot, ddplot, ellipsplot, residualplot,
%                  normqqplot, lsscatter, menureg, menuscoreg,
%                  menupls, menucov, menuls,menuda, menusimca
%

%INTITIALISATION%
if isstruct(out)
    innames=fieldnames(out);
    if strmatch('class',innames,'exact')
        attrib=out.class;
    else
        error('The method had no class identifier.')
    end
    if strmatch('classic',innames,'exact')
        classic=out.classic;
    elseif ismember(attrib,{'CPCA','CPCR','CSIMPLS','MLR','CDA','LS','CSIMCA'})
        classic= 1;
    else
        classic=0;
    end
else
    error('The first inputargument is not a structure.')
end
counter=1;
default=struct('nameplot',0,'labsd',3,'labod',3,'labresd',3,'labmcd',3,'lablts',3,'classic',1);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;
%
if nargin>1
    %
    %placing inputfields in array of strings
    %
    for j=1:nargin-1
        if rem(j,2)~=0
            chklist{i}=varargin{j};
            i=i+1;
        end
    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
    choice=options.nameplot;
    if choice~=0
        ask=0;
    else
        ask=1; %menu of plots
    end
    labsd=options.labsd;
    labod=options.labod;
    labresd=options.labresd;
    labmcd=options.labmcd;
    lablts=options.lablts;
    classicplots=options.classic;
    if ~isstruct(classic) & options.classic==1
        mess=sprintf(['The classical output is not available. Only robust plots will be shown.\n',...
                'Please rerun the preceeding analysis with the option ''classic'' set to 1',...
              '\n if the classical plots are required.']);
        disp(mess)
        classicplots=0;
    end
else
    ask=1; %menu of plots
    choice=0;
    labsd=3;
    labod=3;
    labresd=3;
    labmcd=3;
    lablts=3;
    if ~isstruct(classic) & options.classic==1
       % mess=sprintf(['The classical output is not available. Only robust plots will be shown.\n',...
         %       'Please rerun the preceeding analysis with the option ''classic'' set to 1',...
           %     '\n if the classical plots are required.']);
        %disp(mess)
        classicplots=0;
    else
        classicplots=1;
    end
end

%to initialize the correct menu of plots
exitno=0;
switch attrib
case 'CPCA'
    exitno=4;
    if ismember({char(choice)},{'regdiag','3ddiag','robdist','qqmcd','dd','ellipse','resfit','resindex','qqlts','diag','scatter',...
                'da','simca'})
         error('That kind of plot is not available for this method.')
    end
case 'CPCR'
    exitno=6;
    if ismember({char(choice)},{'robdist','qqmcd','dd','ellipse','resfit','resindex','qqlts','diag','scatter','da','simca'})
         error('That kind of plot is not available for this method.')
    end
case 'LS'
    exitno=7;
     if ismember({char(choice)},{'scree','pcadiag','3ddiag','robdist','qqmcd','dd','ellipse','simca','da'})
         error('That kind of plot is not available for this method.')
    end
case 'MCDREG'
    exitno=3;
     if ismember({char(choice)},{'scree','pcadiag','3ddiag','robdist','qqmcd','dd','ellipse',...
             'resfit','resindex','qqlts','diag','scatter','da','simca'})
         error('That kind of plot is not available for this method.')
    end
case 'MLR'
    exitno=3;
     if ismember({char(choice)},{'scree','pcadiag','3ddiag','robdist','qqmcd','dd','ellipse',...
             'resfit','resindex','qqlts','diag','scatter','da','simca'})
         error('That kind of plot is not available for this method.')
    end
case 'LTS'
    exitno=7;
    if ismember({char(choice)},{'scree','pcadiag','3ddiag',...
            'robdist','qqmcd','dd','ellipse','da','simca'})
         error('That kind of plot is not available for this method.')
    end
    if ~isfield(out,'rd')
        error('Please rerun the LTS-regression again with the option plots equal to 1.')
    end
case {'RAPCA','ROBPCA'}
    exitno=4;
   if ismember({char(choice)},{'regdiag','3ddiag','robdist','qqmcd','dd','ellipse','resfit',...
           'resindex','qqlts','diag','scatter','da','simca'})
        error('That kind of plot is not available for this method.')
   end
case 'RPCR'
      exitno=6;
      if ismember({char(choice)},{'robdist','qqmcd','dd','ellipse','resfit','resindex','qqlts','diag','scatter',...
              'da','simca'})
         error('That kind of plot is not available for this method.')
    end
case {'RSIMPLS','CSIMPLS'}
    exitno=5;
     if ismember({char(choice)},{'scree','robdist','qqmcd','dd','ellipse','resfit','resindex','qqlts','diag','scatter',...
             'da','simca'})
         error('That kind of plot is not available for this method.')
    end
case 'MCDCOV'
    if ismember({char(choice)},{'scree','pcadiag','regdiag','3ddiag','resfit','resindex','qqlts','diag','scatter',...
                'da','simca'})
         error('That kind of plot is not available for this method.')
    end
    exitno=6;
case {'RDA','CDA'}
    exitno=3;
    if ismember({char(choice)},{'scree','pcadiag','3ddiag','robdist','qqmcd','dd','ellipse',...
             'resfit','resindex','qqlts','diag','scatter','regdiag','simca'})
        error('That kind of plot is not available for this method.')
    end
    if size(out.center,2)>2
        disp('Warning: Tolerance ellipses are only drawn for two-dimensional data sets.')
        return
    end
case {'CSIMCA','RSIMCA'}
    exitno=3;
    if ismember({char(choice)},{'scree','pcadiag','3ddiag','robdist','qqmcd','dd','ellipse',...
             'resfit','resindex','qqlts','diag','scatter','regdiag','da'})
        error('That kind of plot is not available for this method.')
    end
    if size(out.pca{1}.P,1) > 3
        disp('Warning: The dimension of the dataset is larger than 3.')
        return
    end
end
if exitno==0
    error(['Your attribute identifier must be one of the following names:',...
            'CPCA, RAPCA, ROBPCA, CPCR, RPCR, LS, LTS, MCDREG, RSIMPLS, CSIMPLS, MCDCOV, CDA,RDA,CSIMCA,RSIMCA'])
end

%plotting what is asked for
if ask==0
    whichplot(out,choice,attrib,exitno,labsd,labod,labresd,labmcd,lablts,classic,classicplots)
else
    %make menu of plots
    while (choice~=exitno)
        switch attrib
            case {'CPCA','ROBPCA','RAPCA'}
                choice=menuscore(out,attrib,exitno,labsd,labod,classic,classicplots);
            case {'MCDREG','MLR'}
                choice=menureg(out,attrib,exitno,labsd,labresd,classic,classicplots);
            case {'COV','MCDCOV'}
                choice=menucov(out,attrib,exitno,labmcd,classic,classicplots);
            case {'RSIMPLS','CSIMPLS'}
                choice=menupls(out,attrib,exitno,labsd,labod,labresd,classic,classicplots);
            case {'CPCR','RPCR'}
                choice=menuscoreg(out,attrib,exitno,labsd,labod,labresd,classic,classicplots);
            case {'LTS','LS'}
                 choice=menuls(out,attrib,exitno,labsd,labresd,lablts,classic,classicplots);
            case{'CDA','RDA'}
                    choice=menuda(out,attrib,exitno,classic,classicplots);
            case{'CSIMCA','RSIMCA'}
                    choice=menusimca(out,attrib,exitno,classic,classicplots);
           end
    end
end

%%%%%%%%%%%%% MAIN FUNCTION %%%%%%%%%%%%%%%
function whichplot(out,plotn,attrib,exitno,labsd,labod,labresd,labmcd,lablts,classic,classicplots)

%Initializing variables

switch attrib

case {'ROBPCA','RAPCA'}
    Xdist=out.sd;
    cutoffX=out.cutoff.sd;
    cutoffO=out.cutoff.od;
    OD=out.od;
    k=out.k;
    L=out.L;
    %classical output
    if isstruct(classic)
        Lcl=classic.L;
        Xdistcl=classic.sd;
        ODcl=classic.od;
        cutoffXcl=classic.cutoff.sd;
        cutoffOcl=classic.cutoff.od;
        attribcl='CPCA';
    end
case 'MCDREG'
    fitted=out.fitted;
    Rdist=out.resd;
    cutoffR=out.cutoff.resd;
    if size(fitted,2)~=1
        multi=1;
    else
        multi=0;
    end
    res=out.res;
    Xdist=out.rd;
    Rsquared=out.rsquared;
    cutoffX=out.cutoff.rd;
    Se=out.cov;
    k=size(out.slope,1);
    %classical output
    if isstruct(classic)
        attribcl=classic.class;
        fittedcl=classic.fitted;
        Rdistcl=classic.resd;
        standrescl=classic.stdres;
        cutoffRcl=classic.cutoff.resd;
        rescl=classic.res;
        Xdistcl=classic.md;
        cutoffXcl=classic.cutoff.md;
        Rsquaredcl=classic.rsquared;
        Secl=classic.cov;
    end
case 'MLR'
    fittedcl=out.fitted;
    Rdistcl=out.resd;
    cutoffRcl=out.cutoff.resd;
    if size(fittedcl,2)~=1
        multi=1;
    else
        multi=0;
    end
    rescl=out.res;
    Xdistcl=out.md;
    cutoffXcl=out.cutoff.md;
    Rsquaredcl=out.rsquared;
    Secl=out.cov;
    k=size(out.slope,1);
    attribcl=out.class;
case 'RPCR'
    fitted=out.fitted;
    Rdist=out.resd;
    cutoffR=out.cutoff.resd;
    if size(fitted,2)~=1
        multi=1;
    else
        multi=0;
    end
    res=out.res;
    Xdist=out.sd;
    Rsquared=out.rsquared;
    cutoffX=out.cutoff.sd;
    cutoffO=out.cutoff.od;
    OD=out.od;
    Se=out.cov;
    k=out.k;
    L=out.robpca.L;
    %classical output
    if isstruct(classic)
        attribcl='CPCR';
        Lcl=classic.cpca.L;
        Xdistcl=classic.sd;
        ODcl=classic.od;
        cutoffOcl=classic.cutoff.od;
        cutoffXcl=classic.cutoff.sd;
        Rdistcl=classic.resd;
        cutoffRcl=classic.cutoff.resd;
    end
case 'LTS'
    resid=out.res;
    fitted=out.fitted;
    scale=out.scale;
    standres=resid/scale;
    n=length(resid);
    if isfield(out,'X')
        x=out.X;
        y=out.y;
    else
        x=0;
        y=0;
    end
    Xdist=out.rd;
    cutoffX=out.cutoff.rd;
    cutoffY=out.cutoff.rd;
    k=size(out.slope,1);
    %classical output
    if isstruct(classic)
        attribcl='LS';
        residcl=classic.res;
        fittedcl=classic.fitted;
        scalecl=classic.scale;
        standrescl=residcl/scalecl;
        Xdistcl=classic.md;
        cutoffXcl=classic.cutoff.md;
        ncl=length(residcl);
        if isfield(classic,'X')
            xcl=classic.X;
            ycl=classic.y;
        else
            xcl=0;
            ycl=0;
        end
    end
case 'LS'
    attribcl=attrib;
    residcl=out.res;
    fittedcl=out.fitted;
    scalecl=out.scale;
    Xdistcl=out.md;
    cutoffXcl=out.cutoff.md;
    standrescl=out.resd;
    ncl=length(residcl);
    if isfield(out,'X')
        xcl=out.X;
        ycl=out.y;
    else
        xcl=0;
        ycl=0;
    end
    k=size(out.slope,1);
case 'CPCA'
    Lcl=out.L;
    Xdistcl=out.sd;
    ODcl=out.od;
    cutoffOcl=out.cutoff.od;
    cutoffXcl=out.cutoff.sd;
    k=out.k;
    attribcl=attrib;
case 'CPCR'
    fitted=out.fitted;
    Lcl=out.cpca.L;
    k=out.k;
    Xdistcl=out.sd;
    ODcl=out.od;
    cutoffOcl=out.cutoff.od;
    cutoffXcl=out.cutoff.sd;
    Rdistcl=out.resd;
    cutoffRcl=out.cutoff.resd;
    if size(fitted,2)~=1
        multi=1; %multivariate analysis
    else
        multi=0; %univariate analysis
        stdrescl=out.resd;
        cutoffRcl=out.cutoff.resd;
    end
    attribcl=attrib;
case 'RSIMPLS'
    fitted=out.fitted;
    Rdist=out.resd;
    cutoffR=out.cutoff.resd;
    if size(fitted,2)~=1
        multi=1;
    else
        multi=0;
    end
    res=out.res;
    Xdist=out.sd;
    cutoffX=out.cutoff.sd;
    cutoffO=out.cutoff.od;
    OD=out.od;
    Se=out.cov;
    k=out.k;
    %classical output
    if isstruct(classic)
        attribcl='CSIMPLS';
        Xdistcl=classic.sd;
        ODcl=classic.od;
        cutoffOcl=classic.cutoff.od;
        cutoffXcl=classic.cutoff.sd;
        Rdistcl=classic.resd;
        cutoffRcl=classic.cutoff.resd;
    end
case 'CSIMPLS'
    fittedcl=out.fitted;
    if size(fittedcl,2)~=1
        multi=1;
    else
        multi=0;
    end
    Xdistcl=out.sd;
    ODcl=out.od;
    k=out.k;
    cutoffXcl=out.cutoff.sd;
    cutoffOcl=out.cutoff.od;
    attribcl=attrib;
    Rdistcl=out.resd;
    cutoffRcl=out.cutoff.resd;
case 'MCDCOV'
    if ~isempty(out.plane)
        disp('Warning (makeplot): The MCD covariance matrix is singular. No plots can be drawn.')
        return
    end
    covar=out.cov;
    md=out.md;
    rd=out.rd;
    if isfield(out,'X')
        data=out.X;
    else
        data=0;
    end
    center=out.center;
    cutoffR=out.cutoff.rd;
    cutoffM=out.cutoff.md;
    if isstruct(classic)
        if isfield(out,'X')
            datacl=out.X;
        else
            datacl=0;
        end
        centercl=out.classic.center;
        covcl=out.classic.cov;
        mdcl=out.classic.md;
        attribcl=out.classic.class;
    end
case 'COV'
   %only available inside the mcdcov function.
case 'CDA'
    if isfield(out,'x')
        xcl=out.x;
        groupcl=out.group;
    else
        xcl=0;
        groupcl=0;
    end
    methodcl=out.method;
    centercl=out.center;
    covcl=out.cov;
    classcl=out.class;
case 'RDA'
    if isfield(out,'x')
        x=out.x; group=out.group;
    else
        x=0;group=0;
    end
    method=out.method;
    center=out.center;
    covar=out.cov;
    class=out.class;
    if isstruct(classic)
        if isfield(out.classic,'x')
            xcl=out.classic.x;
            groupcl=out.classic.group;
        else
            xcl=0;
            groupcl=0;
        end
        methodcl=out.classic.method;
        centercl=out.classic.center;
        covcl=out.classic.cov;
        classcl=out.classic.class;
    end
case 'RSIMCA'
    if isstruct(classic)
        resultcl=out.classic;
    else
        resultcl=0;
    end
end

%determination of the number of plots (robust or/and classic)
if strcmp(attrib,'LS') | strcmp(attrib,'CPCR') | strcmp(attrib,'CPCA') | strcmp(attrib,'MLR') | strcmp(attrib,'CDA') | strcmp(attrib,'CSIMPLS')|strcmp(attrib,'CSIMCA')
    oneplot=1; %only classical plots possible
elseif  (strcmp(attrib,'RAPCA') |strcmp(attrib,'MCDREG') | strcmp(attrib,'RPCR')|strcmp(attrib,'ROBPCA')| strcmp(attrib,'RSIMPLS')| strcmp(attrib,'MCDCOV')|strcmp(attrib,'LTS')|strcmp(attrib,'RDA')|strcmp(attrib,'RSIMCA'))  &  classicplots==0
    oneplot=2; %only robust plots possible
else
    oneplot=0; %both robust and classical plot in the middle of the screen
    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)];
end

switch plotn
case 'all' %all possible plots
    close
    switch attrib
        case {'MCDREG','MLR'}
            if oneplot==0
                figure('Position',pos1)
                regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attrib,labsd,labresd)
                figure('Position',pos2)
                regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,attribcl,labsd,labresd)
            elseif oneplot==1
                figure
                regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,attribcl,labsd,labresd) %multi stond op 2
            elseif oneplot==2
                figure
                regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attrib,labsd,labresd)
            end
        case {'CPCA', 'ROBPCA'}
            if oneplot==0
                figure('Position',pos1)
                screeplot(L,attrib)
                figure('Position',pos2)
                screeplot(Lcl,attribcl)
                pos1(2)=pos1(2)-40;
                pos2(2)=pos2(2)-40;
                figure('Position',pos1)
                scorediagplot(Xdist,OD,k,cutoffX,cutoffO,attrib,labsd,labod)
                figure('Position',pos2)
                scorediagplot(Xdistcl,ODcl,k,cutoffXcl,cutoffOcl,attribcl,labsd,labod)
            elseif oneplot==1
                figure
                screeplot(Lcl,attrib)
                figure
                scorediagplot(Xdistcl,ODcl,k,cutoffXcl,cutoffOcl,attribcl,labsd,labod)
            else
                figure
                screeplot(L,attrib)
                figure
                scorediagplot(Xdist,OD,k,cutoffX,cutoffO,attrib,labsd,labod)
            end
        case {'RSIMPLS','CSIMPLS'}
                if oneplot==0
                    figure('Position',pos1)
                    scorediagplot(Xdist,OD,k,cutoffX,cutoffO,attrib,labsd,labod)
                    figure('Position',pos2)
                    scorediagplot(Xdistcl,ODcl,k,cutoffXcl,cutoffOcl,attribcl,labsd,labod)
                    pos1(2)=pos1(2)-40;
                    pos2(2)=pos2(2)-40;
                    figure('Position',pos1)
                    regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attribcl,labsd,labresd)
                    figure('Position',pos2)
                    regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,attribcl,labsd,labresd)
%3D plot is still highly memory consuming 08/04/04
%                     if exist('OD')
%                         Nihil=(OD <= 1.e-06);
%                         OD(Nihil)=0;
%                         help=OD;
%                     elseif exist('ODcl')
%                         Nihil=(ODcl <= 1.e-06);
%                         ODcl(Nihil)=0;
%                         helpcl=ODcl;
%                     end

%                     pos1(2)=pos1(2)-80;
%                     pos2(2)=pos2(2)-80;
%                     if ~all(help) %in case k=rank(T)=r then OD = zero
%                         figure('Position',pos1)
%                         regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attrib,labsd,labresd)
%                         figure('Position',pos2)
%                         regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,attribcl,labsd,labresd)
%                     else
%                         figure('Position',pos1)
%                         regresdiagplot3d(Xdist,OD,Rdist,cutoffX,cutoffO,cutoffR,k,attrib,multi,labsd,labod,labresd,0)
%                         title(['Robust 3D outlier map based on ',num2str(k),' LV'])
%                         figure('Position',pos2)
%                         regresdiagplot3d(Xdistcl,ODcl,Rdistcl,cutoffXcl,cutoffOcl,cutoffRcl,k,attribcl,multi,labsd,labod,labresd,0)
%                         title(['Classical 3D outlier map based on ',num2str(k),' LV'])
%                     end
                elseif oneplot==1
                    figure
                    scorediagplot(Xdistcl,ODcl,k,cutoffXcl,cutoffOcl,attribcl,labsd,labod)
                    figure
                    regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,attribcl,labsd,labresd)
%3D plot is still highly memory consuming 08/04/04
%                     if exist('ODcl')
%                         Nihil=(ODcl <= 1.e-06);
%                         ODcl(Nihil)=0;
%                         helpcl=ODcl;
%                     end
%                     figure
%                     if ~all(helpcl) %in case k=rank(T)=r then OD = zero
%                         regresdiagplot(Rdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,'CPCR',labsd,labresd)
%                     else
%                         regresdiagplot3d(Xdistcl,ODcl,Rdistcl,cutoffXcl,cutoffOcl,cutoffRcl,k,attribcl,multi,labsd,labod,labresd,0)
%                         title(['Classical 3D outlier map based on ',num2str(k),' LV'])
%                     end
                else
                    figure
                    scorediagplot(Xdist,OD,k,cutoffX,cutoffO,attrib,labsd,labod)
                    figure
                    regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attrib,labsd,labresd)
%3D plot is still highly memory consuming 08/04/04
%                     if exist('OD')
%                         Nihil=(OD <= 1.e-06);
%                         OD(Nihil)=0;
%                         help=OD;
%                     end
%                     figure
%                     if ~all(help) %in case k=rank(T)=r then OD = zero % nog verfijnen is niet volledig nul is zeer klein (<10^-6)
%                         regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attrib,labsd,labresd)
%                     else
%                         regresdiagplot3d(Xdist,OD,Rdist,cutoffX,cutoffO,cutoffR,k,attrib,multi,labsd,labod,labresd,0)
%                         title(['Robust 3D outlier map based on ',num2str(k),' LV'])
%                     end
                end
            case {'CPCR','RPCR'}
            if oneplot==0
                figure('Position',pos1)
                screeplot(L,attrib)
                figure('Position',pos2)
                screeplot(Lcl,attrib)
                pos1(2)=pos1(2)-40;
                pos2(2)=pos2(2)-40;
                figure('Position',pos1)
                scorediagplot(Xdist,OD,k,cutoffX,cutoffO,attrib,labsd,labod)
                figure('Position',pos2)
                scorediagplot(Xdistcl,ODcl,k,cutoffXcl,cutoffOcl,attribcl,labsd,labod)
                pos1(2)=pos1(2)-80;
                pos2(2)=pos2(2)-80;
                figure('Position',pos1)
                regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attribcl,labsd,labresd)
                figure('Position',pos2)
                regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,attribcl,labsd,labresd)
%3D plot is still highly memory consuming 08/04/04
%                 if exist('OD')
%                     Nihil=(OD <= 1.e-06);
%                     OD(Nihil)=0;
%                     help=OD;
%                 elseif exist('ODcl')
%                     Nihil=(ODcl <= 1.e-06);
%                     ODcl(Nihil)=0;
%                     helpcl=ODcl;
%                 end
%                 if ~all(help) %in case k=rank(T)=r then OD = zero
%                     figure('Position',pos1)
%                     regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attrib,labsd,labresd)
%                     figure('Position',pos2)
%                     regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,attribcl,labsd,labresd)
%                 else
%                     figure('Position',pos1)
%                     regresdiagplot3d(Xdist,OD,Rdist,cutoffX,cutoffO,cutoffR,k,attrib,multi,labsd,labod,labresd,0)
%                     title(['Robust 3D outlier map based on ',num2str(k),' LV'])
%                     figure('Position',pos2)
%                     regresdiagplot3d(Xdistcl,ODcl,Rdistcl,cutoffXcl,cutoffOcl,cutoffRcl,k,attribcl,multi,labsd,labod,labresd,0)
%                     title(['Classical 3D outlier map based on ',num2str(k),' LV'])
%                 end
            elseif oneplot==1
                figure
                screeplot(Lcl,attrib)
                figure
                scorediagplot(Xdistcl,ODcl,k,cutoffXcl,cutoffOcl,attrib,labsd,labod)
                figure
                regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,attribcl,labsd,labresd)
%3D plot is still highly memory consuming 08/04/04
%                 if exist('ODcl')
%                     Nihil=(ODcl <= 1.e-06);
%                     ODcl(Nihil)=0;
%                     helpcl=ODcl;
%                 end
%                 figure
%                 if ~all(helpcl) %in case k=rank(T)=r then OD = zero
%                     regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,'CPCR',labsd,labresd)
%                 else
%                     regresdiagplot3d(Xdistcl,ODcl,Rdistcl,cutoffXcl,cutoffOcl,cutoffRcl,k,attribcl,multi,labsd,labod,labresd,0)
%                     title(['Classical 3D outlier map based on ',num2str(k),' LV'])
%                 end
            else
                figure
                screeplot(L,attrib)
                figure
                scorediagplot(Xdist,OD,k,cutoffX,cutoffO,attrib,labsd,labod)
                figure
                regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attrib,labsd,labresd)
%3D plot is still highly memory consuming 08/04/04
%                 figure
%                 if exist('OD')
%                     Nihil=(OD <= 1.e-06);
%                     OD(Nihil)=0;
%                     help=OD;
%                 end
%                 if ~all(help) %in case k=rank(T)=r then OD = zero
%                     regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attrib,labsd,labresd)
%                 else
%                     regresdiagplot3d(Xdist,OD,Rdist,cutoffX,cutoffO,cutoffR,k,attrib,multi,labsd,labod,labresd,0)
%                     title(['Robust 3D outlier map based on ',num2str(k),' LV'])
%                 end
            end
        case 'MCDCOV'
            if oneplot==0
                figure('Position',pos1)
                distplot(rd,cutoffR,attrib,labmcd)
                figure('Position',pos2)
                distplot(md,cutoffM,attribcl,labmcd)
                figure('Position',pos1)
                chiqqplot(rd,size(data,2),attrib)
                figure('Position',pos2)
                chiqqplot(md,size(data,2),attribcl)
                figure
                ddplot(md,rd,cutoffM,attrib,labmcd)
                figure
                if size(data,2)~=2|data==0
                    axes
                    box on
                    text(0.05,0.5,'Tolerance ellipse is only available for bivariate data','color','r')
                else
                    ellipsplot(center,covar,data,rd,labmcd)
                    set(findobj('color','r'),'color','m');
                    hold on
                    ellipsplot(centercl,covcl,data,md,labmcd)
                    set(findobj('color','r'),'color','b','linestyle',':');
                    hold off
                    legend_handles=[findobj('color','m'); findobj('color','b','linestyle',':')];
                    legend(legend_handles,'robust','classical','Location','best')
                end
            elseif oneplot==2
                figure
                distplot(rd,cutoffM,attrib,labmcd)
                figure
                chiqqplot(rd,size(data,2),attrib)
                figure
                ddplot(md,rd,cutoffM,attrib,labmcd)
                figure
                if size(data,2)~=2|data==0
                    axes
                    box on
                    text(0.05,0.5,'Tolerance ellipse is only available for bivariate data','color','r')
                else
                    ellipsplot(center,covar,data,rd,labmcd,'MCDCOV')
                end
            end
        case {'LTS','LS'}
            if oneplot==0
                figure('Position',pos1)
                residualplot(fitted,standres,attrib,'Fitted values','Standardized LTS residual',lablts)
                figure('Position',pos2)
                residualplot(fittedcl,standrescl,attribcl,'Fitted value','Standardized LS residual',lablts)
                pos1(2)=pos1(2)-40;
                pos2(2)=pos2(2)-40;
                figure('Position',pos1)
                residualplot(1:n,standres,attrib,'Index','Standardized LTS residual',lablts)
                figure('Position',pos2)
                residualplot(1:ncl,standrescl,attribcl,'Index','Standardized LS residual',lablts)
                pos1(2)=pos1(2)-80;
                pos2(2)=pos2(2)-80;
                figure('Position',pos1)
                normqqplot(resid,attrib)
                figure('Position',pos2)
                normqqplot(residcl,attribcl)
                pos1(2)=pos1(2)-100;
                pos2(2)=pos2(2)-100;
                figure('Position',pos1)
                if strcmp(attrib,'LTS')
                    regresdiagplot(Xdist,standres,cutoffX,cutoffY,k,0,attrib,labsd,labresd);
                else
                    regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,1,attrib,labsd,labresd);
                end
                figure('Position',pos2)
                if strcmp(attribcl,'LS')
                    regresdiagplot(Xdistcl,standrescl,cutoffXcl,0,k,0,attribcl,labsd,labresd);
                else
                    regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,1,attribcl,labsd,labresd);
                end
                pos1(2)=pos1(2)-120;
                pos2(2)=pos2(2)-120;
                figure('Position',pos1)
                if size(x,2)~=1|x==0
                    axes
                    box on
                    text(0.05,0.5,'Scatter plot is only available for bivariate data','color','r')
                else
                    lsscatter(x,y,fitted,attrib)
                end
                figure('Position',pos2)
                if size(xcl,2)~=1|xcl==0
                    axes
                    box on
                    text(0.05,0.5,'Scatter plot is only available for bivariate data','color','r')
                else
                    lsscatter(xcl,ycl,fittedcl,attribcl)
                end
            elseif oneplot ==1
                figure
                residualplot(fittedcl,standrescl,attribcl,'Fitted value','Standardized LS residual',lablts) ;
                figure
                residualplot(1:ncl,standrescl,attribcl,'Index','Standardized LS residual',lablts) ;
                figure
                normqqplot(residcl,attrib);
                figure
                if strcmp(attrib,'LS')
                    regresdiagplot(Xdistcl,standrescl,cutoffXcl,0,k,0,attrib,labsd,labresd);
                else
                    regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,1,attrib,labsd,labresd);
                end
                figure
                if size(xcl,2)~=1|xcl==0
                    axes
                    box on
                    text(0.05,0.5,'Scatter plot is only available for bivariate data','color','r')
                else
                    lsscatter(xcl,ycl,fittedcl,attribcl)
                end
              elseif oneplot==2
                figure
                residualplot(fitted,standres,attrib, 'Fitted value','Standardized LTS residual',lablts)
                figure
                residualplot(1:n,standres,attrib,'Index','Standardized LTS residual',lablts)
                figure
                normqqplot(resid,attrib)
                figure
                if strcmp(attrib,'LTS')
                    regresdiagplot(Xdist,standres,cutoffX,0,k,0,attrib,labsd,labresd);
                else
                    regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,1,attrib,labsd,labresd);
                end
                figure
                if size(x,2)~=1|x==0
                    axes
                    box on
                    text(0.05,0.5,'Scatter plot is only available for bivariate data','color','r')
                else
                    lsscatter(x,y,fitted,attrib)
                end
            end
        case{'RDA','CDA'}
            if oneplot==0
                    figure('Position',pos1)
                    daplot(x,group,center,covar,class,method)
                    figure('Position',pos2)
                    daplot(xcl,groupcl,centercl,covcl,classcl,methodcl)
                elseif oneplot==1
                    figure
                    daplot(xcl,groupcl,centercl,covcl,classcl,methodcl)
                elseif oneplot==2
                    figure
                    daplot(x,group,center,covar,class,method)
                end
            case{'RSIMCA','CSIMCA'}
                if oneplot==0
                    figure('Position',pos1)
                    simcaplot(result)
                    figure('Position',pos2)
                    simcaplot(resultcl)
                elseif oneplot==1
                    figure
                    simcaplot(resultcl)
                elseif oneplot==2
                    figure
                    simcaplot(result)
                end
        end
    choice=exitno;
case 'scree' %screeplot
    close
    if oneplot==0
        figure('Position',pos1)
        screeplot(L,attrib)
        figure('Position',pos2)
        screeplot(Lcl,attribcl)
    elseif oneplot==1
        screeplot(Lcl,attribcl)
    else
        screeplot(L,attrib)
    end
case 'pcadiag' %diagnostic Scoreplot (Xdist-Odist)
    %close
    if oneplot==0
        figure('Position',pos1)
        scorediagplot(Xdist,OD,k,cutoffX,cutoffO,attrib,labsd,labod)
        figure('Position',pos2)
        scorediagplot(Xdistcl,ODcl,k,cutoffXcl,cutoffOcl,attribcl,labsd,labod)
    elseif oneplot==1
        figure
        scorediagplot(Xdistcl,ODcl,k,cutoffXcl,cutoffOcl,attribcl,labsd,labod)
    else
        figure
        scorediagplot(Xdist,OD,k,cutoffX,cutoffO,attrib,labsd,labod)
    end
case '3ddiag' %3D-plot: highly memory consuming!!!!
    close
    if oneplot==0
        if exist('OD')
            Nihil=(OD <= 1.e-06);
            OD(Nihil)=0;
            help=OD;
        elseif exist('ODcl')
            Nihil=(ODcl <= 1.e-06);
            ODcl(Nihil)=0;
            helpcl=ODcl;
        end
        if ~all(help)%in case k=rank(T)=r then OD = zero
            figure('Position',pos1)
            regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attrib,labsd,labresd)
            figure('Position',pos2)
            regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,attribcl,labsd,labresd)
        else
            figure('Position',pos1)
            regresdiagplot3d(Xdist,OD,Rdist,cutoffX,cutoffO,cutoffR,k,attrib,multi,labsd,labod,labresd,0)
            title(['Robust 3D outlier map based on ',num2str(k),' LV'])
            figure('Position',pos2)
            regresdiagplot3d(Xdistcl,ODcl,Rdistcl,cutoffXcl,cutoffOcl,cutoffRcl,k,attribcl,multi,labsd,labod,labresd,0)
            title(['Classical 3D outlier map based on ',num2str(k),' LV'])
        end
    elseif oneplot==1
        if exist('ODcl')
            Nihil=(ODcl <= 1.e-06);
            ODcl(Nihil)=0;
            helpcl=ODcl;
        end
        if ~all(helpcl) %in case k=rank(T)=r then OD = zero
            regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,attribcl,labsd,labresd)
        else
            regresdiagplot3d(Xdistcl,ODcl,Rdistcl,cutoffXcl,cutoffOcl,cutoffRcl,k,attribcl,multi,labsd,labod,labresd,0)
            title(['Classical 3D outlier map based on ',num2str(k),' LV'])
        end
    else
        if exist('OD')
            Nihil=(OD <= 1.e-06);
            OD(Nihil)=0;
            help=OD;
        end
        if ~all(help) %in case k=rank(T)=r then OD = zero
            regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attrib,labsd,labresd)
        else
            regresdiagplot3d(Xdist,OD,Rdist,cutoffX,cutoffO,cutoffR,k,attrib,multi,labsd,labod,labresd,0)
            title(['Robust 3D outlier map based on ',num2str(k),' LV'])
        end
    end
case {'Idist','robdist'}
    close
    if oneplot==0
        figure('Position',pos1)
        distplot(rd,cutoffM,'MCDCOV',labmcd)
        figure('Position',pos2)
        distplot(md,cutoffM,'COV',labmcd)
        %elseif oneplot==1 not available since 'COV' is only called inside
        %the MCDCOV function
    elseif oneplot==2
        distplot(rd,cutoffM,'MCDCOV',labmcd)
    end
case 'qqmcd'
    close
    if oneplot==0
        figure('Position',pos1)
        chiqqplot(rd,size(data,2),'MCDCOV')
        figure('Position',pos2)
        chiqqplot(md,size(data,2),'COV')
    elseif oneplot==2
        chiqqplot(rd,size(data,2),'MCDCOV')
    end
case 'dd'
    close
    ddplot(md,rd,cutoffM,attrib,labmcd)
case 'ellipse'
    close
    if oneplot==0
        if size(data,2)~=2|data==0
            axes
            box on
            text(0.05,0.5,'Tolerance ellipse plot is only available for bivariate data','color','r')
        else
            ellipsplot(center,covar,data,rd,labmcd)
            set(findobj('color','r'),'color','m');
            hold on
            ellipsplot(centercl,covcl,data,md,labmcd)
            set(findobj('color','r'),'color','b','linestyle',':');
            hold off
            legend_handles=[findobj('color','m'); findobj('color','b','linestyle',':')];
            legend(legend_handles,'robust','classical','Location','best')
        end
    elseif oneplot==2
        if size(data,2)~=2|data==0
            axes
            box on
            text(0.05,0.5,'Tolerance ellipse plot is only available for bivariate data','color','r')
        else
            ellipsplot(center,covar,data,rd,labmcd)
        end
    end
case 'resfit'
    close
    if oneplot==0
        figure('Position',pos1)
        residualplot(fitted,standres,attrib,'Fitted value','Standardized LTS residual',lablts)
        figure('Position',pos2)
        residualplot(fittedcl,standrescl,attribcl,'Fitted value','Standardized LS residual',lablts)
    elseif oneplot==1
        residualplot(fittedcl,standrescl,attribcl,'Fitted value','Standardized LS residual',lablts)
    elseif oneplot==2
        residualplot(fitted,standres,attrib,'Fitted value','Standardized LTS residual',lablts)
    end
case 'resindex'
    close
    if oneplot==0
        figure('Position',pos1)
        residualplot(1:n,standres,attrib,'Index','Standardized LTS residual',lablts)
        figure('Position',pos2)
        residualplot(1:ncl,standrescl,attribcl,'Index','Standardized LS residual',lablts)
    elseif oneplot==1
        residualplot(1:ncl,standrescl,attribcl,'Index','Standardized LS residual',lablts)
    elseif oneplot==2
        residualplot(1:n,standres,attrib,'Index','Standardized LTS residual',lablts)
    end
case 'qqlts'
    close
    if oneplot==0
        figure('Position',pos1)
        normqqplot(resid,'LTS')
        figure('Position',pos2)
        normqqplot(residcl,'LS')
    elseif oneplot==1
        normqqplot(residcl,'LS')
    elseif oneplot==2
        normqqplot(resid,'LTS')
    end
case 'regdiag'
    %close
    if oneplot==0
        figure('Position',pos1)
        if strcmp(attrib,'LTS')
            regresdiagplot(Xdist,standres,cutoffX,0,k,0,attrib,labsd,labresd);
        else
            regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attrib,labsd,labresd);
        end
        figure('Position',pos2)
        if strcmp(attribcl,'LS')
            regresdiagplot(Xdistcl,standrescl,cutoffXcl,0,k,0,attribcl,labsd,labresd);
        else
            regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,attribcl,labsd,labresd);
        end
    elseif oneplot==1
        if strcmp(attrib,'LS')
            regresdiagplot(Xdistcl,standrescl,cutoffXcl,0,k,0,attrib,labsd,labresd);
        else
            regresdiagplot(Xdistcl,Rdistcl,cutoffXcl,cutoffRcl,k,multi,attrib,labsd,labresd);
        end
    elseif oneplot==2
        if strcmp(attrib,'LTS')
            regresdiagplot(Xdist,standres,cutoffX,0,k,0,attrib,labsd,labresd);
        else
            regresdiagplot(Xdist,Rdist,cutoffX,cutoffR,k,multi,attrib,labsd,labresd);
        end
    end
case 'scatter'
    close
    if oneplot==0
        figure('Position',pos1)
        if size(x,2)~=1|x==0
            axes
            box on
            text(0.05,0.5,'Scatter plot is only available for bivariate data','color','r')
        else
            lsscatter(x,y,fitted,attrib)
        end
        figure('Position',pos2)
        if size(xcl,2)~=1|xcl==0
            axes
            box on
            text(0.05,0.5,'Scatter plot is only available for bivariate data','color','r')
        else
            lsscatter(xcl,ycl,fittedcl,attribcl)
        end
    elseif oneplot==1
        if size(xcl,2)~=1|xcl==0
            axes
            box on
            text(0.05,0.5,'Scatter plot is only available for bivariate data','color','r')
        else
            lsscatter(xcl,ycl,fittedcl,attribcl)
        end
    elseif oneplot==2
        if size(x,2)~=1|x==0
            axes
            box on
            text(0.05,0.5,'Scatter plot is only available for bivariate data','color','r')
        else
            lsscatter(x,y,fitted,attrib)
        end
    end
case 'da'
    close
    if oneplot==0
        figure('Position',pos1)
        daplot(x,group,center,covar,class,method)
        figure('Position',pos2)
        daplot(xcl,groupcl,centercl,covcl,classcl,methodcl)
    elseif oneplot==1
        daplot(xcl,groupcl,centercl,covcl,classcl,methodcl)
    elseif oneplot==2
        daplot(x,group,center,covar,class,method)
    end
case 'simca'
    close
    if oneplot==0
        figure('Position',pos1)
        simcaplot(out)
        figure('Position',pos2)
        simcaplot(resultcl)
    elseif oneplot==1
        simcaplot(out)
    elseif oneplot==2
        simcaplot(out)
    end
end

%---------------------------------------------------------------------------------------------------------------------------
function choice= menucov(out,attrib,exitno,labmcd,classic,classicplots)
choice=menu('Covariance plots: ','All','Index plot of the distances','Quantile-Quantile plot of the distances','Distance-distance plot',...
             'Tolerance ellipse (for bivariate data)','Exit');
switch choice
    case 1
        plotn='all';
        whichplot(out,plotn,attrib,exitno,3,3,3,labmcd,3,classic,classicplots)
    case 2
        plotn='Idist';
        whichplot(out,'Idist',attrib,exitno,3,3,3,labmcd,3,classic,classicplots)
    case 3
        plotn='qqmcd';
        whichplot(out,plotn,attrib,exitno,3,3,3,labmcd,3,classic,classicplots)
    case 4
        plotn='dd';
        whichplot(out,plotn,attrib,exitno,3,3,3,labmcd,3,classic,classicplots)
    case 5
        plotn='ellipse';
        whichplot(out,plotn,attrib,exitno,3,3,3,labmcd,3,classic,classicplots)
    case 6
        choice=exitno;
    end
%---------------------------------------------------------------------------------------------------------------------------
function choice = menureg(out,attrib,exitno,labsd,labresd,classic,classicplots)
choice=menu('Regression Plots: ','All ','Regression outlier map', 'Exit ');
switch choice
    case 1
        plotn='all';
        whichplot(out,plotn,attrib,exitno,labsd,3,labresd,3,3,classic,classicplots)
    case 2
        plotn='regdiag';
        whichplot(out,plotn,attrib,exitno,labsd,3,labresd,3,3,classic,classicplots)
    case 3
        choice=exitno;
    end
%---------------------------------------------------------------------------------------------------------------------------
function choice =  menuscore(out,attrib,exitno,labsd,labod,classic,classicplots)
choice=menu('Score Plots: ','All ','Scree plot','Score outlier map', 'Exit ');
switch choice
    case 1
        plotn='all';
        whichplot(out,plotn,attrib,exitno,labsd,labod,3,3,3,classic,classicplots)
    case 2
        plotn='scree';
        whichplot(out,plotn,attrib,exitno,labsd,labod,3,3,3,classic,classicplots)
    case 3
        plotn='pcadiag';
        whichplot(out,plotn,attrib,exitno,labsd,labod,3,3,3,classic,classicplots)
    case 4
        choice=exitno;
    end

%---------------------------------------------------------------------------------------------------------------------------
function choice = menupls(out,attrib,exitno,labsd,labod,labresd,classic,classicplots)
choice=menu('PLS Plots: ','All ','Score outlier map',...
                     'Regression outlier map', '3D outlier map (Highly memory consuming)','Exit ');
switch choice
    case 1
        plotn='all';
        whichplot(out,plotn,attrib,exitno,labsd,labod,labresd,3,3,classic,classicplots)
    case 2
        plotn='pcadiag';
        whichplot(out,plotn,attrib,exitno,labsd,labod,labresd,3,3,classic,classicplots)
    case 3
        plotn='regdiag';
        whichplot(out,plotn,attrib,exitno,labsd,labod,labresd,3,3,classic,classicplots)
    case 4
        plotn='3ddiag';
        whichplot(out,plotn,attrib,exitno,labsd,labod,labresd,3,3,classic,classicplots)
    case 5
        choice=exitno;
    end
%---------------------------------------------------------------------------------------------------------------------------
function choice=menuscoreg(out,attrib,exitno,labsd,labod,labresd,classic,classicplots)
choice=menu('Score and Regression Plots: ','All ','Scree plot','Score outlier map',...
        'Regression outlier map', '3D outlier map  (Highly memory consuming)', 'Exit ');

switch choice
    case 1
        plotn='all';
        whichplot(out,plotn,attrib,exitno,labsd,labod,labresd,3,3,classic,classicplots)
    case 2
        plotn='scree';
        whichplot(out,plotn,attrib,exitno,labsd,labod,labresd,3,3,classic,classicplots)
    case 3
        plotn='pcadiag';
        whichplot(out,plotn,attrib,exitno,labsd,labod,labresd,3,3,classic,classicplots)
    case 4
        plotn='regdiag';
        whichplot(out,plotn,attrib,exitno,labsd,labod,labresd,3,3,classic,classicplots)
    case 5
        plotn='3ddiag';
        whichplot(out,plotn,attrib,exitno,labsd,labod,labresd,3,3,classic,classicplots)
    case 6
        choice=exitno;
    end
%---------------------------------------------------------------------------------------------------------------------------
function choice= menuls(out,attrib,exitno,labsd,labresd,lablts,classic,classicplots)
choice=menu('Residual plots: ', 'All', 'Standardized residuals versus fitted values',...
             'Index plot of standardized residuals','Normal QQplot of residuals',...
             'Diagnostic plot of residuals versus robust distances',...
             'Scatter plot with regression line','Exit');
switch choice
case 1
    plotn='all';
    whichplot(out,plotn,attrib,exitno,labsd,3,labresd,3,lablts,classic,classicplots)
case 2
    plotn='resfit';
    whichplot(out,plotn,attrib,exitno,3,3,3,3,lablts,classic,classicplots)
case 3
    plotn='resindex';
    whichplot(out,plotn,attrib,exitno,3,3,3,3,lablts,classic,classicplots)
case 4
    plotn='qqlts';
    whichplot(out,plotn,attrib,exitno,3,3,3,3,lablts,classic,classicplots)
case 5
    plotn='regdiag';
    whichplot(out,plotn,attrib,exitno,labsd,3,labresd,3,lablts,classic,classicplots)
case 6
    plotn='scatter';
    whichplot(out,plotn,attrib,exitno,3,3,3,3,lablts,classic,classicplots)
case 7
    choice=exitno;
end
%---------------------------------------------------------------------------------------------------------------------------
function choice = menuda(out,attrib,exitno,classic,classicplots)
choice=menu('Discriminant analysis: ','All ','Tolerance ellipses (for bivariate data)', 'Exit ');
switch choice
    case 1
        plotn='all';
        whichplot(out,plotn,attrib,exitno,3,3,3,3,3,classic,classicplots)
    case 2
        plotn='da';
        whichplot(out,plotn,attrib,exitno,3,3,3,3,3,classic,classicplots)
    case 3
        choice=exitno;
end
%------------------------------------------------------------------------------------------------------------------------------
function choice = menusimca(out,attrib,exitno,classic,classicplots)
choice=menu('SIMCA analysis: ','All ','Scatter plot', 'Exit ');
switch choice
    case 1
        plotn='all';
        whichplot(out,plotn,attrib,exitno,3,3,3,3,3,classic,classicplots)
    case 2
        plotn='simca';
        whichplot(out,plotn,attrib,exitno,3,3,3,3,3,classic,classicplots)
    case 3
        choice=exitno;
end

author: marianux

repo: ecg-kit

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


text:
function otherImage = rot90(this, K)
% (Multiple of) 90 deg image rotation; mimicks rot90 in matlab functionality
%
%   Y = MrDataNd()
%   Y.rot90(K)
%
% This is a method of class MrDataNd.
%
% IN
%   K   multiple of 90 degrees used for counterclockwise rotation
%       i.e.    K = 0,1,2 (counterclockwise rotation)
%           or  K = -1,-2,-3,... (clockwise rotation);
%       default : 1
%
% OUT
%   otherImage      MrDataNd where data matrix is rotated and header is updated to
%                   reflect that change
%
% EXAMPLE
%   Y = MrDataNd();
%   Y.rot90(1); % rotate all slices counterclockwise by 90 degrees
%   Y.rot90(-2) % rotate all slices clockwise by 2*90 = 180 degrees
%
%   See also MrDataNd categorical/rot90

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

if nargin < 2
    K = 1;
end

otherImage = this.perform_unary_operation(@(x) rot90(x, K));

% First and second dimension change through rotation around 1, 3 etc.
% multiples of 90 degree...
% doSwapDimensions = mod(K,2) == 1;
% TODO: or shall this be reflected in affineTransformation?
% No, not in affine transformation! We rather take the stance that this is really
% a change of the data and if one wants to adapt the dimInfo, she has to do
% it actively.
doSwapDimensions = false;

if doSwapDimensions
    otherImage.dimInfo.permute([2 1]);
end

author: translationalneuromodeling

repo: tapas


text:
function [vert, face] = read_ply(fn)

% READ_PLY reads triangles, tetraheders or hexaheders from a Stanford *.ply file
%
% Use as
%   [vert, face, prop, face_prop] = read_ply(filename)
%
% Documentation is provided on
%   http://paulbourke.net/dataformats/ply/
%   http://en.wikipedia.org/wiki/PLY_(file_format)
%
% See also WRITE_PLY, WRITE_VTK, READ_VTK

% Copyright (C) 2013-2018, Robert Oostenveld
%
% $Id$

fid = fopen_or_error(fn, 'r');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the file starts with an ascii header
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

line = readline(fid);
if ~strcmp(line, 'ply')
  fclose(fid);
  ft_error('unexpected header line');
end
line = readline(fid);

if ~strncmp(line, 'format', 6)
  fclose(fid);
  ft_error('unexpected header line');
else
  format = line(8:end);
end
line = readline(fid);

if ~strncmp(line, 'element vertex', 14)
  fclose(fid);
  ft_error('unexpected header line');
else
  nvert = str2double(line(16:end));
end
line = readline(fid);

prop = [];
while strncmp(line, 'property', 8)
  tok = tokenize(line);
  prop(end+1).format = tok{2};
  prop(end  ).name   = tok{3};
  line = readline(fid);
end

if ~strncmp(line, 'element face', 12)
  fclose(fid);
  ft_error('unexpected header line');
else
  nface = str2double(line(14:end));
end
line = readline(fid);

if ~strcmp(line, 'property list uchar int vertex_index')    && ...
    ~strcmp(line, 'property list uchar int vertex_indices') && ...
    ~strcmp(line, 'property list uchar uint vertex_index')  && ...
    ~strcmp(line, 'property list uchar uint vertex_indices')
  % the wikipedia documentation specifies vertex_index, but the OPTOCAT files have vertex_indices

  % it would not be very difficult to enhance the reader here with another
  % representation of the faces, i.e. something else than "uchar int"
  fclose(fid);
  ft_error('unexpected header line');
end
line = readline(fid);

while ~strcmp(line, 'end_header')
  line = readline(fid);
end

offset = ftell(fid);
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the file continues with a section of data, which can be ascii or binary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch format

  case 'ascii 1.0'

    fid = fopen_or_error(fn, 'rt');
    fseek(fid, offset, 'cof');
    dat = fscanf(fid,'%f',[numel(prop), nvert])';
    for j=1:length(prop)
      vert.(prop(j).name) = dat(:,j);
    end
    face = zeros(nface,0);
    num = zeros(nface,1);
    for i=1:nface
      % each polygon can have a different number of elements
      num(i) = fscanf(fid,'%f',1);
      face(i,1:num(i)) = fscanf(fid,'%f',num(i));
    end
    fclose(fid);

  case 'binary_little_endian 1.0'

    fid = fopen_or_error(fn, 'rb', 'l');
    fseek(fid, offset, 'cof');

    dat = zeros(nvert,length(prop));
    for i=1:nvert
      for j=1:length(prop)
        dat(i,j) = fread(fid, 1, prop(j).format);
      end % for each property
    end % for each vertex
    for j=1:length(prop)
      switch prop(j).format
        % the format can be one of the following: char uchar short ushort int uint float double int8 uint8 int16 uint16 int32 uint32 float32 float64
        case 'char'
          vert.(prop(j).name) = uint8(dat(:,j));
        case 'uchar'
          vert.(prop(j).name) = uint8(dat(:,j));
        case 'short'
          vert.(prop(j).name) = int16(dat(:,j));
        case 'ushort'
          vert.(prop(j).name) = uint16(dat(:,j));
        otherwise
          vert.(prop(j).name) = dat(:,j);
      end
    end
    clear dat;

    face = zeros(nface,0);
    num  = zeros(nface,1);
    for i=1:nface
      % each polygon can have a different number of elements
      num(i) = fread(fid, 1, 'uint8');
      face(i,1:num(i)) = fread(fid, num(i), 'int32');
    end
    fclose(fid);

  case 'binary_big_endian 1.0'
    % this is exactly the same as the code above, except that the file is opened as big endian

    fid = fopen_or_error(fn, 'rb', 'b');
    fseek(fid, offset, 'cof');

    dat = zeros(nvert,length(prop));
    for i=1:nvert
      for j=1:length(prop)
        dat(i,j) = fread(fid, 1, prop(j).format);
      end % for each property
    end % for each vertex
    for j=1:length(prop)
      switch prop(j).format
        % the format can be one of the following: char uchar short ushort int uint float double int8 uint8 int16 uint16 int32 uint32 float32 float64
        case 'char'
          vert.(prop(j).name) = uint8(dat(:,j));
        case 'uchar'
          vert.(prop(j).name) = uint8(dat(:,j));
        case 'short'
          vert.(prop(j).name) = int16(dat(:,j));
        case 'ushort'
          vert.(prop(j).name) = uint16(dat(:,j));
        otherwise
          vert.(prop(j).name) = dat(:,j);
      end
    end
    clear dat;

    face = zeros(nface,0);
    num  = zeros(nface,1);
    for i=1:nface
      % each polygon can have a different number of elements
      num(i) = fread(fid, 1, 'uint8');
      face(i,1:num(i)) = fread(fid, num(i), 'int32');
    end
    fclose(fid);

  otherwise
    ft_error('unsupported format');
end % switch

% each polygon can have a different number of elements
% mark all invalid entries with nan
if any(num<size(face,2))
  for i=1:nface
    face(i,(num(i)+1):end) = nan;
  end
end % if

if numel(face)>0
  % MATLAB indexes start at 1, inside the file they start at 0
  face = face+1;
end

end % function read_ply

function line = readline(fid)
% read the next line from the ascii header, skip all comment lines
line = fgetl(fid);
while strncmp(line, 'comment', 7)
  line = fgetl(fid);
end
end % function readline

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Z = poly2DApproxLocal( D,  lsX, nrBfsX, lsY, nrBfsY )
%
% Purpose : This function performs local polynomial approximation for
% data which lie on a 2D grid.
%
% This function uses the discrete orthogonal polynomial toolbox DOPbox. The
% toolbox is available at MATLAB file exchange:
%
%   http://www.mathworks.com/matlabcentral/fileexchange/41250
%
% The use of the DOPbox enable the computation with polynomials of very
% high degree. There is no problem to compute a polynomial of degree d =
% 1000, the only disadvantage is the corresponding loss in speed.
%
% Some of the theory associated with this function can be found in the
% publication, Discrete Polynomial Moments and Savitzky-Golay Smoothing,
% available at:
%
% www.waset.org/journals/waset/v48/v48-85.pdf
%
% Use (syntax):
%   Z =  poly2DApproxLocal( D,  lsX, nrBfsX, lsY, nrBfsY )
%
% Input Parameters :
%   D: The 2D data lying on a grid, the function can be applied to any
%   form of data, e.g. surface data or images.
%   nrBfsX: the number of basis functions used in the x direction.
%   lsX: the support length in the x direction.
%   nrBfsY: the number of basis functions used in the y direction.
%   lsY: the support length in the y direction.
%
% Return Parameters :
%   Z: The 2D local polynomial approximation to the data.
%
% Description and algorithms:
%
% References :
%
% Author :  Matther Harker and Paul O'Leary
% Date :    29. Jan 2013
% Version : 1.0
%
% (c) 2013 Matther Harker and Paul O'Leary
% url: www.harkeroleary.org
% email: office@harkeroleary.org
%
% History:
%   Date:           Comment:
%
[ny, nx] = size( D );
%
% generate the rquired polynomial bases
%
xs = (1:nx)';
ys = (1:ny)';
%
Sx = dopApproxLocal( xs, lsX, nrBfsX, 'sparse' );
Sy = dopApproxLocal( ys, lsY, nrBfsY, 'sparse' );
%
%
% Generate the 2D local approximation
%
Z = Sy * D * Sx';

author: Sable

repo: mcbench-benchmarks


text:
% Replace with the location you download the datasets to
data_dir = [getenv('USERPROFILE') '/Dropbox/AAM/test data/'];

% location of datasets, the datasets are acquired from http://ibug.doc.ic.ac.uk/resources/300-W/
dataset_locs = { [data_dir, '/AFW/'];
                 [data_dir, '/helen/trainset/'];
                 [data_dir, '/helen/testset/'];
                 [data_dir, '/ibug/'];
                 [data_dir, 'lfpw/testset/'];
                 [data_dir, '/lfpw/trainset/'];};

bb_locs = {[data_dir, '/Bounding Boxes/bounding_boxes_afw.mat'];
           [data_dir, '/Bounding Boxes/bounding_boxes_helen_trainset.mat'];
           [data_dir, '/Bounding Boxes/bounding_boxes_helen_testset.mat'];
           [data_dir, '/Bounding Boxes/bounding_boxes_ibug.mat'];
           [data_dir, '/Bounding Boxes/bounding_boxes_lfpw_testset.mat'];
           [data_dir, '/Bounding Boxes/bounding_boxes_lfpw_trainset.mat'];};

for d=1:numel(dataset_locs)

    load(bb_locs{d});

    for i=1:numel(bounding_boxes)
        img_name = bounding_boxes{i}.imgName;
        [~, name, ~] = fileparts(img_name);
%         image = imread([dataset_locs{d}, img_name]);
%         imshow(image);
        bbox = bounding_boxes{i}.bb_detector;

        f = fopen([dataset_locs{d}, name, '.txt'], 'w');

        fprintf(f, '%f %f %f %f\r\n', bbox(1), bbox(2), bbox(3), bbox(4));
        fclose(f);

%         rectangle('Position', [bbox(1), bbox(2), bbox(3)-bbox(1), bbox(4)-bbox(2)]);
%         drawnow expose
%         pause(0);
    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function a = reshape(b,varargin)
% Overloaded reshape function for file_array objects
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: reshape.m 1143 2008-02-07 19:33:33Z spm $

if length(struct(b))~=1, error('Can only reshape simple file_array objects.'); end;

args = [];
for i=1:length(varargin),
    args = [args varargin{i}(:)'];
end;
if prod(args)~=prod(b.dim),
    error('To RESHAPE the number of elements must not change.');
end;
a = b;
a.dim = args;


text:
function outStruct = generaldlg(uiStruct,pos,title)
%
% AUTHOR: ???
% PURPOSE: Put up a general dialogue box suited for mrLoadRet
%
% This was uncommented, used in View/Edit Notes call, Annotation
% ... perhaps other places?
% Possibly from rmk?
%
% Probably, we should set the fields and figure to be scalable.
%
% BW, 12.15.00
if isunix
  fontSize = 10;
else
  fontSize = 8;
end
butWidth=10;
topMargin = 1.0;
botMargin = 0.2;
numFields = length(uiStruct);
%scale factors for x and y axis coordinates
xs = 1.45;
ys = 1.6;
yFig = pos(4) * ys;
if exist('title', 'var')
  yStop = zeros(numFields, 1);
  for i=1:numFields, yStop(i)=uiStruct(i).editPos(2) + uiStruct(i).editPos(4); end
  yMax = max(yStop);
  yFig = (yMax + 2) * ys;
end
% Why is resize off?  Because the text boxes don't scale
% when we resize the figure.  There is a way to manage this
% when we build the uistruct.  Sigh.
%
h = figure('MenuBar','none','Units','char',...
    'Position', [pos(1)*xs, pos(2), pos(3)*xs, yFig],...
    'Resize','off',...
    'NumberTitle','off');
bkColor = get(h,'Color');
if exist('title', 'var')
  y = (numFields - topMargin + 5) * ys;
  uicontrol( ...
    'Style', 'text', ...
    'Units', 'char', ...
    'String', title, ...
    'BackgroundColor', bkColor, ...
    'Position', [0, (yMax+0.2)*ys, pos(3)*xs, ys*1.4], ...
    'HorizontalAlignment', 'left', ...
    'FontWeight', 'bold', ...
    'FontSize', 12 ...
    );
end
for uiNum=1:numFields
  %string box
  if ~isempty(uiStruct(uiNum).string)
    stringPos =[uiStruct(uiNum).stringPos(1)*xs,uiStruct(uiNum).stringPos(2)*ys,...
	    uiStruct(uiNum).stringPos(3)*xs,uiStruct(uiNum).stringPos(4)*ys];
    uicontrol('Style','text',...
        'Units','char',...
	'String',uiStruct(uiNum).string,...
	'BackgroundColor',bkColor,...
	'Position',stringPos,...
	'HorizontalAlignment','left',...
	'FontSize',fontSize);
  end
  switch uiStruct(uiNum).style
    case 'edit'
    %edit box
      editPos =[uiStruct(uiNum).editPos(1)*xs,uiStruct(uiNum).editPos(2)*ys,...
	        uiStruct(uiNum).editPos(3)*xs,uiStruct(uiNum).editPos(4)*ys];
       ui_handle(uiNum) = uicontrol('Style','edit',...
     'Units','char',...
	 'String',mat2str(uiStruct(uiNum).value),...
	 'Position',editPos,...
	 'HorizontalAlignment','left',...
	 'FontSize',fontSize);
     case 'text'
      textPos =[uiStruct(uiNum).editPos(1)*xs,uiStruct(uiNum).editPos(2)*ys,...
	        uiStruct(uiNum).editPos(3)*xs,uiStruct(uiNum).editPos(4)*ys];
       ui_handle(uiNum) = uicontrol('Style','text',...
         'Units','char',...
	 'String',mat2str(uiStruct(uiNum).value),...
	 'Position',textPos,...
	 'HorizontalAlignment','left',...
	 'FontSize',fontSize);
    case 'checkbox'
      checkPos =[uiStruct(uiNum).editPos(1)*xs,uiStruct(uiNum).editPos(2)*ys,...
	        uiStruct(uiNum).editPos(3)*xs,uiStruct(uiNum).editPos(4)*ys];
       ui_handle(uiNum) = uicontrol('Style','checkbox',...
         'Units','char',...
	 'String',uiStruct(uiNum).string,...
	 'Value',uiStruct(uiNum).value,...
	 'Position',checkPos,...
	 'HorizontalAlignment','left',...
	 'FontSize',fontSize);
    case 'popupmenu'
      popupPos =[uiStruct(uiNum).editPos(1)*xs,uiStruct(uiNum).editPos(2)*ys,...
	        uiStruct(uiNum).editPos(3)*xs,uiStruct(uiNum).editPos(4)*ys];
      ui_handle(uiNum) = uicontrol('Style','popupmenu',...
         'Units','char',...
	 'String',uiStruct(uiNum).list,...
	 'Value',uiStruct(uiNum).choice,...
	 'Position',popupPos,...
	 'HorizontalAlignment','left',...
	 'FontSize',fontSize);

     otherwise
    disp('Style not recognized');
  end
end
%Display the OK/Cancel buttons
x=1;
y=botMargin;
uicontrol('Style','pushbutton',...
          'String','Cancel',...
	  'Units','char',...
	  'Position',[x*xs,y*ys,butWidth*xs,ys],...
	  'CallBack','uiresume',...
	  'FontSize',fontSize,...
	  'UserData','Cancel');

x = pos(3)-butWidth-1;
uicontrol('Style','pushbutton',...
          'String','OK',...
	  'Units','char',...
	  'Position',[x*xs,y*ys,butWidth*xs,ys],...
	  'CallBack','uiresume',...
	  'FontSize',fontSize,...
	  'UserData','OK');
%let the user select some radio buttons and
%wait for a 'uiresume' callback from OK/Cancel
uiwait
%determine which button was hit.
response = get(gco,'UserData');

%gather the status of the controls if 'OK' was
%selected.  Otherwise return empty matrix.

outStruct = [];
if strcmp(response,'OK')
  for uiNum=1:length(uiStruct)
    switch uiStruct(uiNum).style
      case 'edit'
        %type cast
	if ischar(uiStruct(uiNum).value)
	  val = get(ui_handle(uiNum),'String');
	else
	  val = str2num(get(ui_handle(uiNum),'String'));
	end
      case 'text'
          val = uiStruct(uiNum).value;
      case 'checkbox'
          val = get(ui_handle(uiNum),'Value');
      case 'popupmenu'
          val = char(uiStruct(uiNum).list(get(ui_handle(uiNum),'Value')));
      otherwise
       disp('Style not recognized');
    end
    outStruct = setfield(outStruct,uiStruct(uiNum).fieldName,val);
  end
else
  reply = [];
end
close(h)
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%debug
uiStruct(1).style = 'edit';
uiStruct(1).string = 'Name of your chicken:';
uiStruct(1).stringPos = [1,10,18,1];
uiStruct(1).editPos = [1,9,18,1];
uiStruct(1).value = 'Ralph';
uiStruct(1).fieldName = 'chickenName';

uiStruct(2).style = 'edit';
uiStruct(2).string = 'Name of your cat:';
uiStruct(2).stringPos = [20,10,18,1];
uiStruct(2).editPos = [20,9,18,1];
uiStruct(2).value = 'Fluffy';
uiStruct(2).fieldName = 'catName';

uiStruct(3).style = 'checkbox';
uiStruct(3).string = 'neutered?';
uiStruct(3).stringPos = [1,8,18,1];
uiStruct(3).editPos = [1,7,18,1];
uiStruct(3).value = 1;
uiStruct(3).fieldName = 'catNeuteredFlag';

uiStruct(4).style = 'checkbox';
uiStruct(4).string = 'neutered?';
uiStruct(4).stringPos = [20,8,18,1];
uiStruct(4).editPos = [20,7,18,1];
uiStruct(4).value = 1;
uiStruct(4).fieldName = 'chickenNeuteredFlag';

uiStruct(5).style = 'edit';
uiStruct(5).string = 'IQ:';
uiStruct(5).stringPos = [1,6,9,1];
uiStruct(5).editPos = [10,6,9,1];
uiStruct(5).value = 100.34;
uiStruct(5).fieldName = 'catIQ';

pos = [77,17,40,12];

close all
outStruct = generaldlg(uiStruct,pos);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [res,exact] = AccDot(varargin)
%AccDot       Accurate dot product with faithful, to nearest or K-fold rounding
%
%   res = AccDot(A1,B1,...,An,Bn,K)
%
%res is sum(Ai*Bi)
%   K optional, default K=1 for faithfully rounded result
%   K   =0    rounded to nearest result
%       >1    result in cell array of length K of non-overlapping sequences,
%               sum(res{i})=sum(Ai,Bi) with K-fold accuracy
%       =inf  result in cell array of length K of non-overlapping sequences of
%               sufficient length s.t. sum(res{i})=sum(Ai,Bi) exactly
%       []    interval inclusion of the result; inclusion is best possible if no
%               underflow occurs.
%
%output exact (optional), array of 0/1 of size of result with entry equal to 1
%  iff corresponding dot product exact.
%
%All factors may be real non-interval.
%  The size of every product Ai*Bi must be the same.
%
%One of the factors Ai,Bi may be a cell array. If Ai is a cell array, say,
%  then Ai*Bi is interpreted as sum(Ai{k})*Bi. All entries Ai{k} must be of
%  the same size. This gives a convenient way to compute an approximate inverse
%  of an extremely ill-conditioned matrix only using double precision and our
%  K-fold accurate dot product AccDot. For details and examples with cond(A)~1e65
%  and more, see the INTLAB homepage.
%
%Interval input makes only sense, if all intervals are degenerated (point intervals);
%  therefore omitted.
%
%All factors are first stored into one array to be summed up; summation of
%  components by loops. For matrix input Ai,Bi, there may be interpretation overhead.
%For speed and to fight interpretation overhead, 2 extra matrices are necessary; for
%  the product of to nxn matrices this would be n^3 memory, therefore this is not allowed.
%
%   [res,exact] = AccDot(A1,B1,...,An,Bn,K)
%
%output exact (optional), array of 0/1 of size of result with entry equal to 1
%  iff corresponding dot product exact.
%
%Result is set to NaN if overflow occurs. For simplicity, factors Ai,Bi and
%  products Ai*Bi are limited in absolute value to less or equal 1e300. Small
%  numbers and underflow is treated correctly, rounding to nearest might be off
%  one bit if underflow occurs.
%
%A simple way to compute an inclusion of sum(Ai*Bi) of K-fold accuracy is
%
%   [res,exact] = AccDot(A1,B1,...,An,Bn,K);
%   ResK = intval(res{K});
%   ResK(~exact) = infsup(pred(ResK(~exact)),succ(ResK(~exact)));
%   res{K} = ResK;
%
%This uses faithfully rounded result. Note that the interval is represented by
%   sum(res(1:K-1)) + ResK;
%This interval is generally 2 bits wide unless res(mu,nu)=sum(Ai*Bi)(mu,nu).
%
%A best possible interval including sum(Ai*Bi) is computed by
%
%   res = AccDot(A1,B1,...,An,Bn,[]);
%
%The inclusion is best possible, however, more expensive than the previous approach.
%The gain is at most one bit.
%Maximum number of nonzero elements per sum is limited to 67108862, which
%seems sufficient for Matlab. More elements can be treated, see our paper (Huge length).
%
%Uses various algorithms in
%  S.M. Rump, T. Ogita, S. Oishi: Accurate Floating-point Summation I:
%    Faithful Rounding, SIAM J. Sci. Comput., 31(1):189-224, 2008.
%  S.M. Rump, T. Ogita, S. Oishi: Accurate Floating-point Summation II:
%    Sign, K-fold Faithful and Rounding to Nearest, Siam J. Sci. Comput.,
%    31(2):1269-1302, 2008.
%
%CAUTION: !!! THIS IMPLEMENTATION SUFFERS SEVERELY FROM INTERPRETATION OVERHEAD !!!
%!!! IT IS INCLUDED TO SHOW THE PRINCIPLES OF THE NEW METHODS !!!
%!!! DO NOT USE FOR LARGE DIMENSIONS !!!
%

% written  12/12/05     S.M. Rump
% modified 02/13/06     S.M. Rump  check for complex operands
% modified 08/26/12     S.M. Rump  rounding
% modified 12/23/12     S.M. Rump  typo
%

  e = 1e-30;
  if 1+e==1-e                       % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  % detect parameters
  len = length(varargin);
  if even(len)                          % default accuracy
    K = 1;
  else
    K = varargin{len};
  end

  % check K
  if ~(K>=0) & ~isempty(K)
    error('invalid value K for accuracy in AccDot')
  end

  % detect dimension of result
  Ai = varargin{1};
  if iscell(Ai)
    Ai = Ai{1};
  end
  Bi = varargin{2};
  if iscell(Bi)
    Bi = Bi{1};
  end
  if prod(size(Ai))==1
    [m n] = size(Bi);
  elseif prod(size(Bi))==1
    [m n] = size(Ai);
  else
    m = size(Ai,1);                     % m or n must be 1
    n = size(Bi,2);
  end
  sizeres = [m n];

  % Convert input into A*B+C, where A*B and C is m x n
  % C covers products Ai*Bi where one factor is +/-1, stored as index set indexC
  % index is negative if cofactor is cell array
  A = [];
  B = [];
  indexC = [];                          % summands, cofactor +/- 1
  indexD = [];                          % cell array summands, cofactor +/- 1
  isfull = 1;                           % all factors full

  for i=1:floor(nargin/2)
    % next product Ai*Bi
    Ai = varargin{2*i-1};
    Bi = varargin{2*i};
    if iscell(Ai) & iscell(Bi)
      error('at most factor can be cell array in AccDot')
    end
    if iscell(Ai)
      isfull = isfull & ( ~issparse(Ai{1}) & ~issparse(Bi) );
      sAi = size(Ai{1});
      for k=2:length(Ai)
        if ~isequal(sAi,size(Ai{k}))
          error('sizes in cell array factors in AccDot must coincide')
        end
      end
      sBi = size(Bi);
    elseif iscell(Bi)
      isfull = isfull & ( ~issparse(Ai) & ~issparse(Bi{1}) );
      sAi = size(Ai);
      sBi = size(Bi{1});
      for k=2:length(Bi)
        if ~isequal(sBi,size(Bi{k}))
          error('sizes in cell array factors in AccDot must coincide')
        end
      end
    else
      isfull = isfull & ( ~issparse(Ai) & ~issparse(Bi) );
      sAi = size(Ai);
      sBi = size(Bi);
    end
    % check dimensions
    if prod(sAi)==1
      if ~isequal(sBi,sizeres)
      error('dimensions in AccDot do not match')
      end
    elseif  prod(sBi)==1
      if ~isequal(sAi,sizeres)
        error('dimensions in AccDot do not match')
      end
    elseif sAi(2)~=sBi(1)
      error('dimensions in AccDot do not match')
    end

    % sAi size of factor Ai,  sBi size of factor Bi
    % determine factors
    if iscell(Ai)                         % first factor Ai cell array
      if prod(sBi)==1                     % second factor Bi scalar
        if abs(Bi)==1                     % Bi = +/- 1
          indexD = [ indexD Bi*(2*i-1) ];
        else                              % Bi scalar, but not +/- 1
          if isfull
            Bi_ = sparse(1:n,1:n,Bi);
          else
            Bi_ = Bi*eye(n);
          end
          for k=1:length(Ai)
            A = [ A Ai{k} ];
            B = [ B ; Bi_ ];
          end
        end
      else                                % second factor Bi non-scalar
        if ( prod(sAi)==1 ) & ( ~isequal(sizeres,[1 1]) )
          error('cell array factor cannot be sequence of scalars')
        end
        for k=1:length(Ai)
          A = [ A Ai{k} ];
          B = [ B ; Bi ];
        end
      end
    elseif iscell(Bi)                     % second factor Bi cell array
      if prod(sAi)==1                     % first factor Ai scalar
        if abs(Ai)==1                     % Ai = +/- 1
          indexD = [ indexD Ai*(2*i) ];
        else                              % Ai scalar, but not +/- 1
          if isfull
            Ai_ = sparse(1:m,1:m,Ai);
          else
            Ai_ = Ai*eye(m);
          end
          for k=1:length(Bi)
            A = [ A Ai_ ];
            B = [ B ; Bi{k} ];
          end
        end
      else                                % first factor Ai non-scalar
        if ( prod(sBi)==1 ) & ( ~isequal(sizeres,[1 1]) )
          error('cell array factor cannot be sequence of scalars')
        end
        for k=1:length(Bi)
          A = [ A Ai ];
          B = [ B ; Bi{k} ];
        end
      end
    else                                  % none of Ai,Bi cell array
      if prod(sAi)==1                     % first factor Ai scalar
        if abs(Ai)==1                     % Ai = +/- 1
          indexC = [ indexC Ai*(2*i) ];
        else
          if isfull
            A = [ A Ai*eye(m) ];
          else
            A = [ A sparse(1:m,1:m,Ai) ]; % Ai scalar, but not +/- 1
          end
          B = [ B ; Bi ];
        end
      elseif prod(sBi)==1                 % second factor Bi scalar
        if abs(Bi)==1                     % Bi = +/- 1
          indexC = [ indexC Bi*(2*i-1) ];
        else
          A = [ A Ai ];                   % Bi scalar, but not +/- 1
          if isfull
            B = [ B ; sparse(1:n,1:n,Bi) ];
          else
            B = [ B ; Bi*eye(n) ];
          end
        end
      else                                % both factors Ai,Bi non-scalar
        A = [ A Ai ];
        B = [ B ; Bi ];
      end
    end
  end

  % check interval input
  if isa(A,'intval') | isa(B,'intval')
    error('input of AccDot must not be of type intval')
  end

  % check for complex input
  if ( ~isreal(A) ) | ( ~isreal(B) )
    error('input of AccDot must be real')
  end

  % check huge input
  if any(any(abs(A)>1e300)) | any(any(abs(B)>1e300))
    error('for simplicity all entries and intermediate products limited to 1e300')
  end
  if ~isempty(indexC)
    for k=1:length(indexC)
      if any(any(varargin{abs(indexC(k))}>1e300))
        error('for simplicity all entries and intermediate products limited to 1e300')
      end
    end
  end
  kD = 0;
  if ~isempty(indexD)
    for k=1:length(indexD)
      DD = varargin{abs(indexD(k))};
      kD = kD + length(DD);
      for kk=1:length(DD)
        if any(any(DD{kk}>1e300))
          error('for simplicity all entries and intermediate products limited to 1e300')
        end
      end
    end
  end

  % transpose A:  result A.'*B
  A = A.';

  factor = 134217729;               % splitting factor 2^27+1
  % [A1,A2] = Split(A)
  C = factor*A;
  A1 = C - ( C - A );               % upper part of A
  A2 = A - A1;                      % A = A1+A2 exact splitting
  % [B1,B2] = Split(B)
  C = factor*B;
  B1 = C - ( C - B );               % upper part of B
  B2 = B - B1;                      % B = B1+B2 exact splitting

  % initialize result
  if isempty(K)
    rescell = 0;
  else
    rescell = ( K>1 );
  end
  if isfull
    if rescell
      res{1} = zeros(m,n);
    else
      res = zeros(m,n);
    end
  else
    if rescell
      res{1} = sparse([],[],[],m,n);
    else
      res = sparse([],[],[],m,n);
    end
  end
  if ( K>1 ) & ~isinf(K)
    for k=2:K
      res{k} = res{1};
    end
  end
  if isempty(K)
    res = intval(res);
  end
  if nargout==2
    if isfull
      exact = zeros(m,n);
    else
      exact = sparse([],[],[],m,n);
    end
  end

  % compute result
  sA = size(A,1);
  sC = length(indexC);
  sD = length(indexD);
  summands = ( sC+sD~=0 );
  % underflow constants
  phi = 2^643;                                % factor sqrt(eps^-4 eta^-1) for Ai,Bi
  Fbound = 2^(-968);                          % lower bound eps^-2 eta for factors
  for i=1:m
    for j=1:n
      if sA~=0
        x1 = A1(:,i);
        x2 = A2(:,i);
        y1 = B1(:,j);
        y2 = B2(:,j);
        AiBip = A(:,i).*B(:,j);
        r = x2.*y2 - ((( AiBip - x1.*y1 ) - x2.*y1 ) - x1.*y2 );
        if any(abs(AiBip)>1e300)              % AiBip is vector
          error('for simplicity all entries and intermediate products limited to 1e300')
        end
      else
        AiBip = [];
        r = [];
      end
      if summands
        Ci = zeros(sC+kD,1);
        for k=1:sC
          if indexC(k)>0
            Ci(k) = varargin{indexC(k)}(i,j);
          else
            Ci(k) = -varargin{-indexC(k)}(i,j);
          end
        end
        kkk = 0;
        for k=1:sD
          DD = varargin{abs(indexD(k))};
          if indexD(k)>0
            for kk=1:length(DD)
              kkk = kkk+1;
              Ci(sC+kkk) = DD{kk}(i,j);
            end
          else
            for kk=1:length(DD)
              kkk = kkk+1;
              Ci(sC+kkk) = -DD{kk}(i,j);
            end
          end
        end
      else
        Ci = [];
      end
      p = [ AiBip ; Ci ];
      % execute one ExtractVector
      np = nnz(p);                          % initialization
      mu = full(max(abs(p)));               % abs(p_i) <= mu; full: avoid matlab bug
      if ( np==0 ) | ( mu==0 )              % no or only zero summands
        R = 0;                              % result exactly zero in any rounding
        p = [];
      else
        Ms = 2^nextpow2(np+2);              % np+2 <= 2^M
        if Ms^2*eps>1
          error('vector length n too large for AccDot; Huge n not implemented')
        end
        sigma = Ms*2^nextpow2(mu);          % first extraction unit
        if isinf(sigma) | isnan(sigma)      % overflow, could be avoided with scaling
          if iscell(res)
            res{1}(i,j) = NaN;
          else
            res(i,j) = NaN;
          end
          exact(i,j) = 0;
          continue
        else
          q = ( sigma + p ) - sigma;        % [R,p] = ExtractVector(sigma,p);
          R = sum(q);                       % sum of leading terms
          p = p - q;                        % remaining terms
        end
      end
      if K<=1                               % nearest of faithful rounding
        [res(i,j),ext,R,p] = AccSum( [ p ; r ] , K , 0 , R );
        if res(i,j)==realmin                % take care of underflow, original sum = R + sum(p)
          index = find( abs(AiBip)<Fbound );      % small parts of AiBip
          p = [ p ; -AiBip(index) ; -r(index) ];  % subtract small parts
          [res(i,j),ext] = ...
            AccSum( scaleinput(A(index,i),B(index,j),p,phi) , K , 0 , phi*(phi*R) );
          if ( nargout==2 ) & ext
            exact(i,j) = diam( (intval(res(i,j))/phi) / phi == 0 );
          end
          res(i,j) = (res(i,j)/phi) / phi;  % avoid overflow
        end
      elseif isempty(K)                       % interval result
        [rr,ext,R,p] = AccSum( [ p ; r ] , 0 , -1 , R );
        if rr==realmin                        % take care of underflow, original sum = R + sum(p)
          index = find( abs(AiBip)<Fbound );  % small parts of AiBip
          p = [ p ; -AiBip(index) ; -r(index) ];  % subtract small parts
          [rr,ext] = AccSum( scaleinput(A(index,i),B(index,j),p,phi) , 0 , -1 , phi*(phi*R) );
          if ext
            res(i,j) = intval(rr);
          else
            res(i,j) = infsup(rr,succ(rr));
          end
          res(i,j) = ( res(i,j)/phi ) / phi;
        else
          res(i,j) = infsup(rr,succ(rr));
        end
      else                                          % K-fold accuracy
        p = [ p ; r ];
        scale = 0;
        cont = 1;
        k = 0;
        ext(i,j) = 0;
        while ( k<K ) & cont
          k = k+1;
          [rr,ext,R,p] = AccSum( p , 1 , 0 , R );
          if ~scale & ( rr==realmin )               % take care of underflow, original sum = R + sum(p)
            scale = 1;                              % only one scaling
            index = find( abs(AiBip)<Fbound );      % small parts of AiBip
            p = [ p ; -AiBip(index) ; -r(index) ];  % subtract small parts
            [rr,ext,R,p] = ...                      % scaled result
              AccSum( scaleinput(A(index,i),B(index,j),p,phi) , 1 , 0 , phi*(phi*R) );
          end
          rrr = rr;
          if scale
            rr = ( rr/phi )/phi;
          end
          if isinf(K) & ( k>length(res) ) & ( rr~=0 )
            if isfull
              res{k} = zeros(m,n);
            else
              res{k} = sparse([],[],[],m,n);
            end
          end
          if rr~=0
            res{k}(i,j) = rr;
          end
          if abs(rr)<=realmin                       % result exact
            cont = 0;
          end
          if isnan(rr)
            res{k}(i,j) = NaN;
            cont = 0;
          end
          if ( nargout==2 ) & ext
            exact(i,j) = ( diam( (intval(rrr)/phi) / phi )==0 );
          end
        end
      end
    end
  end

  if rndold
    setround(rndold)
  end

function [res,exact,R,p] = AccSum(p,K,rnd,rho)
%Simplified version of AccSum :
%  rnd only valid for K==0
%  no overflow can occur
%  indices for possible underflow stored
%  vector length not too large
%  stops if sigma too small
%  ~(K>1)
%

  res = 0;
  exact = 1;
  R = 0;
  if isempty(p)
    return
  end

  % input real, compute sum(p)
  % rounding to nearest
  if K==0                             % rnd maybe -1 or 0 or +1
    % let S:=sum(p)
    [res,dummy,R,p] = AccSum(p,1,0,rho);   % S = res + R + sum(p)  for new p
    if res==realmin                   % take care of underflow
      return
    end
    [delta,dummy,R,p] = AccSum(p,1,0,R);   % delta + R + sum(p) = S - res  for new R,p
    if delta==realmin                 % scale; otherwise sign(delta) = sign(S-res)
      return
    end
    if delta==0                       % result exact: res=S in any rounding
      return
    end
    exact = 0;                        % result not exact

    if rnd~=0                         % rounding downwards or upwards
      if sign(delta)==rnd             % res on wrong side
        if rnd==1
          res = succ(res);
        else
          res = pred(res);
        end
      end
      return
    end

    % sign(delta) = sign(S-res)
    % compute nearest neighbor resp in direction sign(delta), rnd=0
    if delta>0
      resp = succ(res);
    else
      resp = pred(res);
    end
    mu = (resp-res)/2;
    if abs(delta)>abs(mu)
      res = resp;
    elseif abs(delta)==abs(mu)
      delta = AccSum(p,1,0,R);
      if delta==0
        res = res + mu;
      elseif sign(delta)==sign(mu)
        res = resp;
      end
    end
    return
  end

  % the standard case: real vector input, K=1
  if issparse(p)
    n = nnz(p);                         % initialization
  else
    n = length(p);
  end
  tau1 = 0;
  tau2 = 0;
  mu = full(max(abs(p)));               % abs(p_i) <= mu
  if ( n==0 ) | ( mu==0 )               % no or only zero summands
    res = rho;                          % result exactly zero in any rounding
    return
  end
  Ms = 2^nextpow2(n+2);                 % n+2 <= 2^M
  sigma = Ms*2^nextpow2(mu);            % first extraction unit
  phi = 2^(-53)*Ms;                     % factor to decrease sigma
  factor = 2*phi*Ms;                    % factor for sigma check

  % underflow constant
  sigmamin = 2^(-968);                  % lower bound eps^-2 eta for sigma
  t = rho;
  while 1
    if sigma<=sigmamin                  % sigma too small, but t+tau exact
      res = realmin;                    % indicates underflow
      R = t;                            % original sum = R + sum(p)
      return
    end
    q = ( sigma + p ) - sigma;          % [tau,p] = ExtractVector(sigma,p);
    tau = sum(q);                       % sum of leading terms
    p = p - q;                          % remaining terms
    tau1 = t + tau;                     % new approximation
    if abs(tau1)>=factor*sigma
      tau2 = tau - ( tau1 - t );        % [tau1,tau2] = FastTwoSum(t,tau)
      res = tau1 + ( tau2 + sum(p) );   % faithfully rounded final result
      R = tau2 - ( res - tau1 );        % only for K-fold result
      if nargout>=2
        exact = ( R==0 ) & ~any(p);
      end
      return
    end
    t = tau1;                           % sum t+tau exact
    if t==0                             % accelerate case sum(p)=0
      [res,exact,R,p] = AccSum(p(p~=0),K,0,0);  % recursive call, zeros eliminated
      return
    end
    sigma = phi*sigma;                  % new extraction unit
  end

function p = scaleinput(Ai,Bi,Ci,phi)
% scale  p = phi*Ai * phi*Bi  for column vector input

  factor = 134217729;               % splitting factor 2^27+1
  % [A1,A2] = Split(phi*Ai)
  Ai = phi*Ai;
  C = factor*Ai;
  A1 = C - ( C - Ai );            % upper part of Ai
  A2 = Ai - A1;                   % Ai = A1+A2 exact splitting
  % [B1,B2] = Split(phi*Bi)
  Bi = phi*Bi;
  C = factor*Bi;
  B1 = C - ( C - Bi );            % upper part of Bi
  B2 = Bi - B1;                   % Bi = B1+B2 exact splitting
  % phi^2*A(:,i)*B(:,j) = p+r, no underflow possible
  p = Ai.*Bi;
  p = [ p ; A2.*B2 - ((( p - A1.*B1 ) - A2.*B1 ) - A1.*B2 ) ; phi*(phi*Ci) ];

repo: OpenMAS