function result=csimpls(x,y,varargin)

%CSIMPLS performs Partial Least Squares regression using the SIMPLS 
% algorithm of de Jong (1993).
%
% Reference: 
%    de Jong, S. (1993),
%    "SIMPLS: an alternative approach to Partial Least Squares Regression", 
%    Chemometrics and Intelligent Laboratory Systems, 18, 251-263.
%
% Required input arguments:
%            x : Data matrix of the explanatory variables
%                (n observations in rows, p variables in columns)
%            y : Data matrix of the response variables
%                (n observations in rows, q variables in columns)
%     
% Optional input arguments:
%            k : Number of components to be used.  
%                (default = min([9,rank([x,y]),floor(n/2),p])). 
%        plots : If equal to one, a menu is shown which allows to draw several plots,
%                such as a score outlier map and a regression outlier map. (default)
%                If 'plots' is equal to zero, all plots are suppressed.
%                See also makeplot.m
%
% I/O: result=csimpls(x,y,'k',10);
%
% The output of CSIMPLS is a structure containing:
%
%   result.slope     : Classical slope estimate
%   result.int       : Classical intercept estimate
%   result.fitted    : Classical prediction vector
%   result.res       : Classical residuals
%   result.cov       : Estimated variance-covariance matrix of the residuals 
%   result.M         : Classical center of the matrix [X;Y]
%   result.T         : Classical scores
%   result.weights.r : Classical simpls weights
%   result.weights.p : Classical simpls weights
%   result.Tcov      : Classical covariance matrix of the scores
%   result.k         : Number of components used in the regression
%   result.sd        : Classical score distances
%   result.od        : Classical orthogonal distances
%   result.resd        : Residual distances (when there are several response variables).
%                      If univariate regression is performed, it contains the standardized residuals.
%   result.cutoff    : Cutoff values for the score, orthogonal and residual distances.
%   result.flag      : The observations whose orthogonal distance is larger than result.cutoff.od
%                      (orthogonal outliers => result.flag.od) and/or whose residual distance is 
%                      larger than result.cutoff.resd (bad leverage points/vertical outliers => result.flag.resd)
%                      can be considered as outliers and receive a flag equal to zero (=> result.flag.all).
%                      The regular observations, including the good leverage points, receive a flag 1.
%   result.class     : 'CSIMPLS'
%
% 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/04/2003 
%Last revision: 20/05/2008

[n,p1]=size(x);
[n2,q1]=size(y);
if n~=n2
    error('The response variables and the predictor variables have a different number of observations.')
end
kmin=min([9,rank(x),floor(n/2),p1]);
default=struct('plots',1,'k',kmin);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;   
counter=1;
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

[xcentr,cx]=mcenter(x);
[ycentr,cy]=mcenter(y);
out.M=[cx cy];%[xcentr ycentr];
sigmaxy=xcentr'*ycentr;
for i=1:options.k
    sigmayx=sigmaxy';
    if q1>p1        
        [RR,LL]=eig(sigmaxy*sigmayx);  
        [LL,I]=greatsort(diag(LL));
        rr=RR(:,I(1));
        qq=sigmayx*rr; 
        qq=qq/norm(qq);
    else
        if q1==1
            qq = 1;
            rr = sigmaxy;
        else
            [QQ,LL]=eig(sigmayx*sigmaxy);
            [LL,I]=greatsort(diag(LL));
            qq=QQ(:,I(1));
            rr=sigmaxy*qq;
        end
    end
    tt=xcentr*rr;
    nttc=norm(tt);
    rr=rr/nttc;
    tt=tt/nttc;
    qq=ycentr'*tt;
    uu=ycentr*qq;
    pp=xcentr'*tt;
    vv=pp;
    if i>1 												
        vv = vv -v*(v'*pp);
    end
    vv = vv/norm(vv);
    sigmaxy = sigmaxy - vv*(vv'*sigmaxy);
    v(:,i)=vv;
    q(:,i)=qq;
    t(:,i)=tt;
    u(:,i)=uu;
    p(:,i)=pp;
    r(:,i)=rr;
end

%Second Stage : Classical Regression and transformation

b=r*q';
int=cy-cx*b;

%classical output
out.T=t;
out.weights.p=p;
out.weights.r=r;
out.coef=[b; int];
out.b=b;
out.int=int;
out.yhat=x*b+repmat(int,n,1);
out.res=y-out.yhat;
out.covar=cov(out.res);
if q1==1
    out.stdres=out.res./sqrt(out.covar);
end
out.k=options.k;
STTm=sum((y-repmat(mean(y),length(y),1)).^2);
SSE=sum(out.res.^2);
out.rsquare=1-SSE/STTm;
out.class='CSIMPLS';

%calculating classical distances in x-space
out.Tcov=cov(t);
out.sd=sqrt(libra_mahalanobis(t,zeros(1,out.k),'cov',out.Tcov))';
out.cutoff.sd=sqrt(chi2inv(0.975,out.k));

%calculating classical orthogonal distances
xt=t*p';
tempo=xcentr-xt;
for i=1:n
    out.od(i,1)=norm(tempo(i,:));
end
r=rank(x);
if out.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

%calculating residual distances
if q1>1
    if (-log(det(out.covar)/(p1+q1-1)))>50
        out.resd='singularity';
    else
       cen=zeros(q1,1);
       out.resd=sqrt(libra_mahalanobis(out.res,cen','cov',out.covar))';
   end
else    %here q==1
    out.resd=out.stdres; %standardized residuals
end
out.cutoff.resd=sqrt(chi2inv(0.975,q1));

%Computing flags
out.flag.od=out.od<=out.cutoff.od;
out.flag.resd=abs(out.resd)<=out.cutoff.resd;
out.flag.all=(out.flag.od & out.flag.resd);


result=struct('slope',{out.b}, 'int',{out.int},'fitted',{out.yhat},'res',{out.res}, 'cov',{out.covar},...
    'M', {out.M},'T',{out.T} ,'weights', {out.weights},'Tcov',{out.Tcov},'k',{out.k},'sd', {out.sd},'od',{out.od},'resd',{out.resd},...
    '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 
end