function result = csimca(x,group,varargin);

%CSIMCA performs the SIMCA method. This is a classification
% method on a data matrix x with a known group structure. On each group a 
% robust PCA analysis is performed. Afterwards a classification
% rule is developped to determine the assignment of new observations. 
%
% Required input arguments:
%          x : training data set (matrix of size n by p).
%      group : column vector containing the group numbers of the training
%              set x. For the group numbers, any strict positive integer is
%              allowed.
%
% Optional input arguments:
%              k : Is a vector with size equal to the number of groups, or otherwise 0. Represents the number
%                  of components to be retained in each group. (default = 0.)
%         method : Indicates which classification rule is wanted. `1' results in rule (R1)
%                  based on the scaled orthogonal and score distances. `2' corresponds with
%                  (R2) based on the squared scaled orthogonal and score distances. Default is 2. 
%          gamma : Represents the value(s) used in the classification rule: weight gamma is given to the od's,
%                  weight (1-gamma) to the sd's. (default = 0.5).
%     misclassif : String which indicates how to estimate the probability of
%                  misclassification. It can be based on the training data ('training'), 
%                  a validation set ('valid'), or cross-validation ('cv'). Default is 'training'.
% membershipprob : Vector which contains the membership probability of each
%                  group (sorted by increasing group number). These values are used to
%                  obtain the total misclassification percentage. 
%          valid : If misclassif was set to 'valid', this field should contain 
%                  the validation set (a matrix of size m by p).
%     groupvalid : If misclassif was set to 'valid', this field should contain the group numbers
%                  of the validation set (a column vector).
%     predictset : Contains a new data set (a matrix of size mp by p) from which the 
%                  class memberships are unknown and should be predicted.  
%          plots : If equal to 1, one figure is created with the training data and the
%                  boundaries for each group. This plot is
%                  only available for trivariate (or smaller) data sets. For technical reasons, a maximum 
%                  of 6 groups is allowed. Default is one.
%       plotspca : If equal to one, a score diagnostic plot is
%                  drawn (default). If 'plots' is equal to zero, this plot is suppressed.
%                  See also makeplot.m
%          labsd : The 'labsd' observations with largest score distance are
%                  labeled on the diagnostic plot. (default = 3)
%          labod : The 'labod' observations with largest orthogonal distance are
%                  labeled on the diagnostic plot. default = 3) 
%
% Options for advanced users (input comes from the program RSIMCA.m with option 'classic' = 1):
%
%   weightstrain : The weights for the training data. Corresponds to the flagtrain from RSIMCA. (default = 1) 
%   weightsvalid : The weights for the validation data. Corresponds to the flagvalid from RSIMCA. (default = 1)
%
% I/O: result=csimca(x,group,'method',1,'misclassif','training',...
%                  'membershipprob',proportions,'valid',y,'groupvalid',groupy,'plots',0);
%
%  The user should only give the input arguments that have to change their default value.
%  The name of the input arguments needs to be followed by their value.
%  The order of the input arguments is of no importance.
%
% Examples: out=csimca(x,group,'method','1')
%           out=csimca(x,group,'plots',0)
%           out=csimca(x,group,'valid',y,'groupvalid',groupy)
%
% The output is a structure containing the following fields:
%        result.assignedgroup : If there is a validation set, this vector contains the assigned group numbers
%                               for the observations of the validation set. Otherwise it contains the
%                               assigned group numbers of the original observations based on the discriminant rules.
%                result.pca   : A cell containing the results of the different PCA analysis on the training sets.
%               result.method : String containing the method used to obtain
%                               the discriminant rules (either 1 for 'R1' or 2 for 'R2'). This
%                               corresponds to the input argument method. 
%            result.flagtrain : Observations from the training set whose score distance and/or orthogonal distance
%                               exceeds a certain cut-off value can be considered as outliers and receive a flag equal 
%                               to zero. The regular observations receive a flag 1. (See also robpca.m)
%            result.flagvalid : Observations from the validation set whose score distance and/or orthogonal distance 
%                               exceeds a certain cut-off value can be considered as outliers and receive a
%                               flag equal to zero. The regular observations receive a flag 1. 
%                               If there is no validation set, this field is equal to zero.
%         result.grouppredict : If there is a prediction set, this vector contains the assigned group numbers
%                               for the observations of the prediction set. 
%          result.flagpredict : Observations from the new data set (predict) whose robust distance (to the center of their group)
%                               exceeds a certain cut-off value can be considered as overall outliers and receive a
%                               flag equal to zero. The regular observations receive a flag 1. 
%                               If there is no prediction set, this field is
%                               equal to zero.
%       result.membershipprob : A vector with the membership probabilities.  If no priors are given, they are estimated 
%                               as the proportions of observations in the training set.
%           result.misclassif : String containing the method used to estimate the misclassification probabilities
%                               (same as the input argument misclassif)
%     result.groupmisclasprob : A vector containing the misclassification probabilities for each group.
%       result.avemisclasprob : Overall probability of misclassification (weighted average of the misclassification
%                               probabilities over all groups).
%                result.class : 'CSIMCA'
%                    result.x : The training data set (same as the input argument x) (only in output when p<=3).
%                result.group : The group numbers of the training set (same as the input argument group) (only in output when p<=3).
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at: 
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Karlien Vanden Branden  
% Last Update: 05/07/2005
%

if nargin<2
    error('There are too few input arguments.')
end

% assigning default-values
[n,p]=size(x);
if size(group,1)~=1
    group=group';
end
if n ~= length(group)
    error('The number of observations is not the same as the length of the group vector!')
end
g=group;
counts=tabulate(g); %contingency table (outputmatrix with 3 colums): value - number - percentage 
[lev,levi,levj]=unique(g);
if ~all(counts(:,2)) %some groups have zero values, omit those groups
    disp(['Warning: group(s) ', num2str(counts(counts(:,2)==0,1)'), 'are empty']);
    empty=counts(counts(:,2)==0,:);
    counts=counts(counts(:,2)~=0,:);
else
    empty=[];
end
ng=size(counts,1);
proportions = zeros(ng,1);
y=0; %initial values of the validation data set and its groupsvector
groupy=0;
labsd = 3;
labod = 3;
counter=1;
gamma = 0.5;
k = zeros(ng,1);
weightstrain = ones(1,n);
weightsvalid = 0;
default=struct('k',k,'method',2,'gamma',0.5,'misclassif','training',...
    'membershipprob',proportions,'valid',y,'groupvalid',groupy,'plots',1,'plotspca',0,'labsd',labsd,...
    'labod',labod,'weightstrain',weightstrain,'weightsvalid',0,'predictset',[]);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;
%reading the user's input 
if nargin>2
    %
    %placing inputfields in array of strings
    %
    for j=1:nargin-2
        if rem(j,2)~=0
            chklist{i}=varargin{j};
            i=i+1;
        end
    end 
    %
    %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
end

%Checking gamma
gamma = options.gamma;
if gamma >1 | gamma <0
    error('An inappropriate number for gamma is given. A correct value lies between 0 and 1.');
end

%Checking prior (>0 )
prior=options.membershipprob;
if size(prior,1)~=1
    prior=prior';
end
epsilon=10^-4;
if sum(prior) ~=0 & (any(prior < 0) | (abs(sum(prior)-1)) > epsilon)
   error('Invalid membership probabilities.')
end
if length(prior)~=ng
    error('The number of membership probabilities is not the same as the number of groups.')
end

%%%%%%%%%%%%%%%%%%MAIN FUNCTION %%%%%%%%%%%%%%%%%%%%%
%Checking if a validation set is given
if strmatch(options.misclassif, 'valid','exact') 
    if options.valid==0
        error(['The misclassification error will be estimated through a validation set',...
                'but no validation set is given!'])
    else
        validx = options.valid;
        validgroup = options.groupvalid;
        if size(validx,1)~=length(validgroup)
            error('The number of observations in the validation set is not the same as the length of its group vector!')
        end
        if size(validgroup,1)~=1
            validgroup = validgroup';
        end
        countsvalid=tabulate(validgroup);
        countsvalid=countsvalid(countsvalid(:,2)~=0,:);
        if size(countsvalid,1)==1
            error('The validation set must contain observations from more than one group!')
        elseif any(ismember(empty,countsvalid(:,1)))
            error(['Group(s) ' ,num2str(empty(ismember(empty,countsvalid(:,1)))), 'was/were empty in the original dataset.'])
        end
    end
    if (length(options.weightsvalid) == 1) | (length(options.weightsvalid)~=size(validx,1))
        options.weightsvalid = ones(size(validx,1),1);
    end
elseif options.valid~=0
    validx = options.valid;
    validgroup = options.groupvalid;
    if size(validx,1) ~= length(validgroup)
        error('The number of observations in the validation set is not the same as the length of its group vector!')
    end
    if size(validgroup,1)~=1
        validgroup = validgroup';
    end
    options.misclassif='valid';
    countsvalid=tabulate(validgroup);
    countsvalid=countsvalid(countsvalid(:,2)~=0,:);
    if size(countsvalid,1)==1
        error('The validation set must contain more than one group!')
    elseif any(ismember(empty,countsvalid(:,1)))
        error(['Group(s) ' , num2str(empty(ismember(empty,countsvalid(:,1)))), ' was/were empty in the original dataset.'])
    end
    if (length(options.weightsvalid) == 1) | (length(options.weightsvalid)~=size(validx,1))
        options.weightsvalid = ones(size(validx,1),1);
    end
end

model.counts = counts(:,2);
model.x = x;
model.group = group;

labsd = floor(max(0,min(options.labsd,n)));
labod = floor(max(0,min(options.labod,n)));

%CSIMCA:
%PRINCIPAL COMPONENT ANALYSIS
%Perform PCA on each group separately:
%   a) if k is not given: decide on the optimal number of components using CV.
%   b) if k is given: perform cpca with the optimal number of components.

for iClass = 1:ng
    indexgroup = find(g==iClass);
    groupi = x(indexgroup,:); 
    if options.k(iClass) == 0
        disp(['A scree plot is drawn for group ',num2str(iClass),'.'])
    end
    model.result{iClass} = cpca(groupi,'k',options.k(iClass),'plots',options.plotspca,...
        'labsd',labsd,'labod',labod);
    model.flag(1,indexgroup) = model.result{iClass}.flag.all;
end

%CLASSIFICATION
%Discriminant rule based on the training set x
[odsc,sdsc] = testmodel(model,x);
finalgrouptrain = assigngroup(odsc,sdsc,options.method,ng,gamma); 
    
if sum(prior) == 0
    result1.prior = counts(:,3)'./100;
else
    result1.prior = prior;
end

%Compute scaled orthogonal and scaled score distances for the validation set
if strmatch(options.misclassif,'valid','exact') 
    [odsc,sdsc] = testmodel(model,validx);
    finalgroup = assigngroup(odsc,sdsc,options.method,ng,gamma);
elseif strmatch(options.misclassif,'cv','exact')  %use cv
    model.k = k;
    [odsc,sdsc] = leave1out(model);
    finalgroup = assigngroup(odsc,sdsc,options.method,ng,gamma);
end

switch options.misclassif
case 'valid'
    [v,vi,vj]=unique(validgroup);
    odscgroup = [];
    sdscgroup = [];
    for iClass = 1:ng
        indexgroup = find(validgroup == iClass);
        odscgroup = [odscgroup;odsc(indexgroup,iClass)];
        sdscgroup = [sdscgroup;sdsc(indexgroup,iClass)];
    end 
    weightsvalid=zeros(length(odscgroup),1); 
    weightsvalid(((odscgroup <= 1) & (sdscgroup <= 1)))=1;
    for igamma = 1:length(gamma)
        for iClass=1:ng 
            misclas(iClass)=sum((validgroup(options.weightsvalid==1)==finalgroup(options.weightsvalid==1,igamma)') & ...
                (validgroup(options.weightsvalid==1)==repmat(v(iClass),1,sum(options.weightsvalid))));
            ingroup(iClass) = sum((validgroup(options.weightsvalid == 1) == repmat(v(iClass),1,sum(options.weightsvalid))));
        end
        misclas = (1 - (misclas./ingroup));
        misclasprobpergroup(igamma,:)=misclas;
        misclas=misclas.*result1.prior;
        misclasprob(igamma)=sum(misclas);    
    end
case 'training'
    for igamma = 1:length(gamma)
        for iClass = 1:ng
            result1.misclas(iClass) = sum((group(options.weightstrain==1)==finalgrouptrain(options.weightstrain==1,igamma)')& ...
                (group(options.weightstrain==1)==repmat(lev(iClass),1,sum(options.weightstrain))));
            result1.ingroup(iClass) = sum((group(options.weightstrain == 1) == repmat(lev(iClass),1,sum(options.weightstrain))));
        end
        misclas = (1 - (result1.misclas./result1.ingroup));
        misclasprobpergroup(igamma,:) = misclas;
        misclas = misclas.*result1.prior;
        misclasprob(igamma) = sum(misclas);
    end
    weightsvalid=0;%only available with validation set
    finalgroup = finalgrouptrain;
case 'cv' 
    for igamma = 1:length(gamma)
        for iClass=1:ng 
            misclas(iClass)=sum((group(options.weightstrain==1)==finalgroup(options.weightstrain==1,igamma)') & ...
                (group(options.weightstrain==1)==repmat(lev(iClass),1,sum(options.weightstrain))));
            ingroup(iClass) = sum((group(options.weightstrain == 1) == repmat(lev(iClass),1,sum(options.weightstrain))));
        end
        misclas = (1 - (misclas./ingroup));
        misclasprobpergroup(igamma,:)=misclas;
        misclas=misclas.*result1.prior;
        misclasprob(igamma)=sum(misclas);    
    end
    weightsvalid=0; %only available with validation set
end

if ~isempty(options.predictset)
    [odscpredict,sdscpredict] = testmodel(model,options.predictset);
    finalgrouppredict = assigngroup(odscpredict,sdscpredict,options.method,ng,gamma)';
    weightspredict = max((odscpredict <= 1) & (sdscpredict <= 1),[],2)';
else
    finalgrouppredict = 0;
    weightspredict = 0;
end   

%Output structure
result = struct('assignedgroup',{finalgroup'},'pca',{model.result},'method',options.method,...
    'flagtrain',{model.flag},'flagvalid',{weightsvalid'},'grouppredict',finalgrouppredict,'flagpredict',weightspredict,...
    'membershipprob',{result1.prior},'misclassif',{options.misclassif},'groupmisclasprob',{misclasprobpergroup},...
    'avemisclasprob',{misclasprob},'class',{'CSIMCA'},'x',x,'group',group);

if size(x,2)>3
    result=rmfield(result,{'x','group'});
end

%Plots:
try
    if options.plots
        makeplot(result)
    end
catch %output must be given even if plots are interrupted 
    %> delete(gcf) to get rid of the menu 
end

%---------------------------
%Leave-One-Out procedure 

function [odsc,sdsc] = leave1out(model)

nClass = length(model.result);

 for iClass = 1:nClass
    index = 1;
    indexgroup = find(model.group == iClass);
    teller_if_lus = 0;
    groupi = model.x(indexgroup,:);
    for i = 1:model.counts(iClass)
        groupia = removal(groupi,i,0);
        GRes = model.result;
        GRes{iClass} = cpca(groupia,'k',model.result{iClass}.k,'plots',0);
        %Calculate for each class the sd and the od for the observation that was left out:
        for jClass = 1:nClass
            dataicentered = model.x(indexgroup(index),:)-GRes{jClass}.M;
            scorei = dataicentered*GRes{jClass}.P;
            dataitilde = scorei*GRes{jClass}.P';
            sd(indexgroup(index),jClass) = sqrt(scorei*(diag(1./GRes{jClass}.L))*scorei');
            od(indexgroup(index),jClass) = norm(dataicentered-dataitilde);
            if GRes{jClass}.cutoff.od ~= 0  
                odsc(indexgroup(index),jClass) = od(indexgroup(index),jClass)/GRes{jClass}.cutoff.od;
            else
                odsc(indexgroup(index),jClass) = 0;
            end
            sdsc(indexgroup(index),jClass) = sd(indexgroup(index),jClass)/GRes{jClass}.cutoff.sd;
        end
        index = index + 1;
    end
end

%--------------------------------
function [odsc,sdsc] = testmodel(model,validx);

%Apply the given model on the test data to obtain different 
%orthogonal distances and score distances.

nClass = length(model.result);
n = size(validx,1);

for jClass = 1:nClass
    for index = 1:n
        out{jClass} = model.result{jClass};
        dataicentered = validx(index,:)-out{jClass}.M;
        scorei = dataicentered*out{jClass}.P;
        dataitilde = scorei*out{jClass}.P';
        sd(index,jClass) = sqrt(scorei*(diag(1./out{jClass}.L))*scorei');
        od(index,jClass) = norm(dataicentered-dataitilde);
        if out{jClass}.cutoff.od ~= 0  
            odsc(index,jClass) = od(index,jClass)/out{jClass}.cutoff.od;
        else
            odsc(index,jClass) = 0;
        end
        sdsc(index,jClass) = sd(index,jClass)/out{jClass}.cutoff.sd; 
    end
end

%-------------------------
function result = assigngroup(odsc,sdsc,method,nClass,gamma);

%Obtain the group assignments for given od's and sd's.

if method == 1 
    sd = sdsc;
    od = odsc;
elseif method == 2 
    sd = sdsc.^2;
    od = odsc.^2;
end

for igamma = 1:length(gamma)
    tdist = gamma(igamma).*od + (1-gamma(igamma)).*sd;
    for i = 1:size(od,1)
        result(i,igamma) = find(tdist(i,:) == min(tdist(i,:)));
    end
end