function result=cpca(x,varargin)

%CPCA performs a classical principal components analysis 
% on the data set x. The loadings correspond with the eigenvectors of the
% classical covariance matrix of x.
% 
% Required input argument:
%      x : Data matrix (observations in the rows, variables in the
%          columns)
%
% Optional input argument: 
%      k : Number of principal components to compute. If k is missing, 
%          a screeplot is plotted allowing the selection of
%          the number of principal components.
%  plots : If equal to one (default), a menu is shown which allows to draw several plots,
%          such as a score outlier map. 
%          If 'plots' is equal to zero, all plots are suppressed.
%
% I/O: result=cpca(x,'k',2)
%
% The output is a structure containing 
%
%    result.P       : Classical loadings (eigenvectors of covariance matrix)
%    result.L       : Classical eigenvalues (eigenvalues of covariance matrix)
%    result.M       : Classical mean of the columns of x
%    result.T       : Classical scores 
%    result.k       : Number of (chosen) principal components 
%    result.sd      : Score (mahalanobis) distances within the classical PCA subspace
%    result.od      : Orthogonal distances to the classical PCA subspace 
%    result.cutoff  : Cutoff values for the score and orthogonal distances
%    result.flag     : The observations whose score distance is larger than result.cutoff.sd (==> result.flag.sd)
%                      or whose orthogonal distance is larger than result.cutoff.od (==> result.flag.od)
%                      can be considered as outliers and receive a flag equal to zero (result.flag.all).
%                      The regular observations receive a flag 1.
%    result.class   : 'CPCA' 
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at: 
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Sabine Verboven 
% Last Update: 05/04/2003

default=struct('plots',1,'k',0);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;   
counter=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
end
[n,p]=size(x);
% using the optimal (least time consuming) algorithm
if n<p
	[P,T,L,r,Xc,out.M]=kernelEVD(x);
else
	[P,T,L,r,Xc,out.M]=classSVD(x);
end

if options.k==0
    screeplot(L,'CPCA');
    k=input('How many principal components would you like to retain? ');
else
    k=options.k;
end
out.k=k;
out.P=P(:,1:out.k);
out.T=T(:,1:out.k);
out.L=L(1:out.k)';

% Mahalanobis distance in classical PCA subspace
Tclas=Xc*out.P;
out.sd=sqrt(libra_mahalanobis(Tclas,zeros(size(Tclas,2),1),'cov',out.L'));
out.cutoff.sd=sqrt(chi2inv(0.975,out.k));

% Orthogonal distances to classical PCA subspace
Xtilde=Tclas*out.P';
Cdiff=Xc-Xtilde;
for i=1:n
    out.od(i)=norm(Cdiff(i,:));
end
r=rank(x);
if k~=r
    m=mean(out.od.^(2/3));
    s=sqrt(var(out.od.^(2/3)));
    out.cutoff.od = sqrt(norminv(0.975,m,s)^3); 
else
    out.cutoff.od=0;
end

%Computing flags
if k~=p
    out.flag.od=(out.od<=out.cutoff.od)';
    out.flag.sd=(out.sd<=out.cutoff.sd)';
    out.flag.all=(out.flag.od)&(out.flag.sd);
else
    out.flag.od=(out.od<=out.cutoff.od)';
    out.flag.sd=(out.sd<=out.cutoff.sd)';
    out.flag.all=out.flag.sd;
end
out.class='CPCA';

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

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