function result=rsimpls(x,y,varargin)

%RSIMPLS is a 'Robust method for Partial Least Squares Regression based on the
% SIMPLS algorithm'. It can be applied to both low and high-dimensional predictor variables x
% and to one or multiple response variables y. It is resistant to outliers in the data.
% The RSIMPLS algorithm is built on two main stages. First, a matrix of scores is derived 
% based on a robust covariance criterion (see robpca.m),
% and secondly a robust regression is performed based on the results from ROBPCA. 
% 
% The RSIMPLS method is described in: 
%    Hubert, M., and Vanden Branden, K. (2003),
%    "Robust Methods for Partial Least Squares Regression",
%    Journal of Chemometrics, 17, 537-549.
%
% To select the number of components in the regression model, a robust RMSECV (root mean squared
% error of cross validation) curve is drawn, based on a fast algorithm for
% cross-validation. This approach is described in:
%
%    Engelen, S., Hubert, M. (2005),
%    "Fast model selection for robust calibration methods",
%    Analytica Chimica Acta, 544, 219-228.
%
% 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(rank([x,y]),kmax)). If k is not specified,
%                it can be selected using the option 'rmsecv'. 
%         kmax : Maximal number of components to be used (default = 9).
%                If k is provided, kmax does not need to be specified, unless k is larger
%                than 9.                 
%        alpha : (1-alpha) measures the fraction of outliers the algorithm should 
%                resist. Any value between 0.5 and 1 may be specified. (default = 0.75)
%            h : (n-h+1) measures the number of observations the algorithm should 
%                resist. Any value between n/2 and n may be specified. (default = 0.75*n)
%                Alpha and h may not both be specified.
%       rmsecv : If equal to zero and k is not specified, a robust R-squared curve
%                is plotted and an optimal k value can be chosen. (default)
%                If equal to one and k is not specified, a robust component selection-curve is plotted and
%                an optimal k value can be chosen. This curve computes a combination of
%                the cross-validated root mean squared error and the robust residual 
%                sum of squares for k=1 to kmax. Rmsecv and k may not both
%                be specified.
%        rmsep : If equal to one, the robust RMSEP-value (root mean squared error of
%                prediction) for the model with k components (default = 0). This value
%                is automatically given if rmsecv = 1.
%        plots : If equal to one, a menu is shown which allows to draw several plots,
%                such as a robust score outlier map and a regression
%                outlier map are drawn. (default)
%                If the input argument 'classic' is equal to one, the classical
%                diagnostic plots are drawn as well.
%                If 'plots' is equal to zero, all plots are suppressed.
%                See also makeplot.m
%        labsd : The 'labsd' observations with largest score distance are
%                labeled on the outlier map (default = 3)
%        labod : The 'labod' observations with largest orthogonal distance are
%                labeled on the outlier map (default = 3)    
%      labresd : The 'labresd' observations with largest residual distance are
%                labeled on the outlier map (default = 3)   
%      classic : If equal to one, the classical SIMPLS analysis will be performed as well
%                (see also csimpls.m). (default = 0)
%
% Options for advanced users:
%           kr : Total number of components used by the ROBPCA method.
%                We advise to use kr=k+q. (default) 
%        kmaxr : Maximal number of components used by the ROBPCA method. 
%                default = min(kmax+q,rank([x,y]))
%  plotsrobpca : If equal to one, a robust score outlier map from ROBPCA is drawn. 
%                If the input argument 'classic' is equal to one, the classical
%                outlier map is drawn as well.
%                If 'plotsrobpca' is equal to zero, all plots are suppressed. (default)
%           st : Indicates the current stage of the algorithm for cross-validation (RMSECV)  
%                default = 0 -> performs all the stages of the algorithm and computes all the 
%                parameters for the full model. 
%                st=1, robpca is still performed, but when 
%                st=2, the algorithm proceeds based on the previous knownledge of the output from robpca. 
%          out : Is an empty structure, but is constructed while performing the cross-validation.
%                (see rrmse.m)
%
% I/O: result=rsimpls(x,y,'k',k,'kmax',10,'alpha',0.75,'h',h,'rmsecv',0,'rmsep',0,...
%       'plots',1,'labsd',3,'labod',3,'labresd',3,'classic',1,'kr',kr,...
%       'kmaxr',kmaxr,'plotsrobpca',0,'st',0,'out',[]);
%  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:
%   rsimpls(x,y,'k',5,'plots',1);
%   rsimpls(x,y,'classic',1,'rmsecv',1);
%
% The output of RSIMPLS is a structure containing:
%
%   result.slope     : Robust slope estimate
%   result.int       : Robust intercept estimate
%   result.fitted    : Robust fitted values
%   result.res       : Robust residuals
%   result.cov       : Estimated variance-covariance matrix of the residuals 
%   result.T         : Robust scores
%   result.weights.r : Robust simpls weights
%   result.weights.p : Robust simpls weights
%   result.Tcenter   : Robust center of the scores
%   result.Tcov      : Robust covariance matrix of the scores 
%   result.rsquared  : Robust R-squared value for the optimal k
%   result.rcs       : Robust Component Selection Criterion:
%                      This is a matrix with kmax columns. The first row contains the approximate
%                      R-squared values (for k=1,...,kmax) and the second row the square root of the weighted
%                      residuals sum of squares. This is equal to the RCS-value with 
%                      gamma = 0 (see Engelen and Hubert (2005) for the definition). 
%                      If the input argument rmsecv = 1, the third row contains the RCS-value for
%                      gamma = 0.5 and the fourth for gamma = 1. The last one is equal to the robust 
%                      cross-validated RMSE values. 
%                      Note that all the entries in this matrix depend on the choice of kmax. 
%   result.rmsep     : Robust RMSEP value 
%   result.k         : Number of components used in the regression
%   result.h         : The quantile h used throughout the algorithm
%   result.sd        : Robust score distances
%   result.od        : Robust 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 (result.cutoff.sd), orthogonal 
%                      (result.cutoff.od) and residual distances (result.cutoff.resd).
%                      We use 0.975 quantiles of the chi-squared distribution.
%   result.flag      : The observations whose score distance is larger than 
%                      'result.cutoff.sd' receive a flag 'result.flag.sd' equal
%                      to zero (good leverage points). Otherwise 'result.flag.sd'
%                      is equal to one. 
%                      The components 'result.flag.od' and 'result.flag.resd' are
%                      defined analogously, and determine the orthogonal outliers, 
%                      resp. the bad leverage points/vertical outliers. 
%                      The observations with 'result.flag.od' and 'result.flag.resd'
%                      equal to zero, can be considered as calibration outliers and receive
%                      'result.flag.all' equal to zero. The regular observations and the good leverage
%                      points have 'result.flag.all' equal to one.
%   result.class     : 'RSIMPLS'
%   result.classic   : If the inputargument 'classic' is equal to 1, this structure
%                      contains results of the classical SIMPLS analysis. (see also csimpls.m)
%   results.robpca   : The results of robpca on [X,Y]. 
%
% 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
% Version date: 07/04/2004 
% Last update: 04/08/2006

%
%initialization with defaults
%
if rem(nargin,2)~=0
    error('Number of input arguments must be even!');
end
[n,p1]=size(x);
[n2,q1]=size(y);
z=[x,y];
rx=rank(x);
rz=rank(z);
if n~=n2
    error('The response variables and the predictor variables have a different number of observations.')
end
niter=100;counter=1;mcd=0;alfa=0.75;
kmax=min([9,rx,floor(n/2),p1]);
kmaxr=min([kmax+q1,rz]);
h=floor(2*floor((n+kmaxr+1)/2)-n+2*(n-floor((n+kmaxr+1)/2))*alfa);
labsd=3;labod=3;labresd=3;
plotsrobpca=0;plots=1;k=0;
kr=k+q1;
st=0;rmsecv=0;
out=[];classic=0;rmsep=0;rmsep_value=nan;rmsecv_value = nan;rsquared_value = nan;rss_value = nan;
default=struct('alpha',alfa,'h',h,'labsd',labsd,'labod',labod,'labresd',labresd,'k',k,'kr',kr,...
    'plotsrobpca',plotsrobpca,'plots',plots,'kmax',kmax,'st',st,...
    'out',out,'rmsecv',rmsecv,'classic',classic,'rmsep',rmsep,...
    'kmaxr',kmaxr,'rmsep_value',rmsep_value,'rmsecv_value',rmsecv_value,'rsquared_value',...
    rsquared_value,'rss_value',rss_value);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;
%
if nargin>2
    %
    %placing inputfields in array of strings
    %
    for j=1:nargin-3
        if rem(j,2)~=0
            chklist{i}=varargin{j};
            i=i+1;
        end
    end 
    dummy=sum(strcmp(chklist,'h')+2*strcmp(chklist,'alpha'));
    switch dummy
        case 0 %Take on default values  
            options.alpha=alfa;%0.75;
            options.h=h;
        case 3
            error('Both input arguments alpha and h are provided. Only one is required.')
    end
    %
    %Checking which default parameters have to be changed
    % and keep them in the structure 'options'.
    %
    while counter<=IN 
        index=strmatch(list(counter,:),chklist,'exact');%contains the users input one by one
        if ~isempty(index) %in case of similarity
            for j=1:nargin-3 %searching the index of the accompanying field
                if rem(j,2)~=0 %fieldnames are placed on odd index
                    if strcmp(chklist{index},varargin{j})
                        I=j;
                    end
                end
            end
            options=setfield(options,chklist{index},varargin{I+1});
            index=[];
        end
        counter=counter+1;
    end
    options.h=floor(options.h);
    options.kmax=floor(options.kmax);
    options.k=floor(options.k);
    options.kmaxr=floor(options.kmaxr);
    options.kr=floor(options.kr);
    kmax=min([options.kmax,floor(n/2),rz,p1]);
    kmaxr=max([min([options.kmaxr,rz]),kmax+q1]);
    k=min(options.k,kmax);
    while k<0
        k=input(['The number of components can not be negative.\n'...
        'How many principal components would you like to retain?\n']);
    end
    if any(strcmp(chklist,'kr'))
        kr=floor(max([options.kr,k+q1]));
    else
        kr=k+q1;
    end
    if dummy==1 %checking inputvariable h
        if options.h-floor(options.h)~=0
            mess=sprintf('Attention (rsimpls.m): h must be an integer. \n');
            disp(mess)
        end
        if kr==0
            if options.h<floor((n+kmaxr+1)/2) 
                options.h=floor((n+kmaxr+1)/2);
                mess=sprintf(['Attention (rsimpls.m): h should be larger than (n+kmaxr+1)/2.\n',...
                        'It is set to its minimum value',num2str(options.h)]);
                disp(mess)
            end
            options.alpha=options.h/n;
        else
            if options.h<floor((n+kr+1)/2) 
                options.h=floor((n+kr+1)/2);
                mess=sprintf(['Attention (rsimpls.m): h should be larger than (n+kr+1)/2.\n',...
                        'It is set to its minimum value',num2str(options.h)]);
                disp(mess)
            end
            options.alpha=options.h/n;
        end
        if options.h>n
            options.alpha=0.75;
            if kr==0
                options.h=floor(2*floor((n+kmaxr+1)/2)-n+2*(n-floor((n+kmaxr+1)/2))*options.alpha);
            else
                options.h=floor(2*floor((n+kr+1)/2)-n+2*(n-floor((n+kr+1)/2))*options.alpha);
            end    
            mess=sprintf(['Attention (rsimpls.m): h should be smaller than n. \n',...
                    'It is set to its default value ',num2str(options.h)]);
            disp(mess)
        end
    elseif dummy==2
        if options.alpha < 0.5
            options.alpha=0.5;
            mess=sprintf(['Attention (rsimpls.m): Alpha should be larger than 0.5. \n',...
                    'It is set to 0.5.']);
            disp(mess)
        end
        if options.alpha > 1
            options.alpha=0.75;
            mess=sprintf(['Attention (rsimpls.m): Alpha should be smaller than 1.\n',...
                    'It is set to 0.75.']);
            disp(mess)
        end
        if kr==0
            options.h=floor(2*floor((n+kmaxr+1)/2)-n+2*(n-floor((n+kmaxr+1)/2))*options.alpha);
        else
            options.h=floor(2*floor((n+kr+1)/2)-n+2*(n-floor((n+kr+1)/2))*options.alpha);
        end  
    end
    h=options.h;alfa=options.alpha;labsd=max(0,min(floor(options.labsd),n));
    dummyh = strcmp(chklist,'h');
    dummykmax = strcmp(chklist,'kmax');
    if all(dummyh == 0) && any(dummykmax)
        h = floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*alfa);
    end
    labod=max(0,min(floor(options.labod),n));labresd=max(0,min(floor(options.labresd),n));
    plotsrobpca=options.plotsrobpca;plots=options.plots;
    st=options.st;
    out=options.out;
    rmsecv=options.rmsecv;
    classic=options.classic;
    rmsep=options.rmsep;
    rmsep_value = options.rmsep_value;
    rmsecv_value = options.rmsecv_value;
    rmsecv = options.rmsecv;
    rsquared_value = options.rsquared_value;
    rss_value = options.rss_value;
end

if q1==1 && k>=(h-2)
    mess=sprintf(['Attention (rsimpls.m): The number of components, k = ',num2str(k),...
            '\n is larger than our recommended maximum value of k = ',num2str(h-2)-1,'.']);
    disp(mess)
elseif q1>1 && k>=((h/q1)-(q1/2)-0.5)
    mess=sprintf(['Attention (rsimpls.m): The number of components, k = ',num2str(k),...
            '\n is larger than our recommended maximum value of k = ',num2str(floor((h/q1)-(q1/2)-1.5)),'.']);
    disp(mess)
end

%
%MAIN PART
%
% selection of number of components
if k == 0
    if rmsecv == 0
        [R2,final]=rsquared(x,y,kmax,'RSIMPLS',options.h);
        rss = final.rss;
        k = final.k;
        result=rsimpls(x,y,'k',k,'kr',k+q1,'h',h,'rmsecv',0,'rmsep',options.rmsep,'plots',0,'classic',classic,'rsquared_value',R2,'rss_value',rss); 
    else
        out=rrmse(x,y,h,kmax,'RSIMPLS',1);
        R2 = out.R2;
        rss = out.rss;
        k=out.k;
        rmsecv_value = out.rmsecv;
        pred = rrmse(x,y,h,kmax,'RSIMPLS',0,k,out.weight,out.res);
        result=rsimpls(x,y,'k',k,'kr',k+q1,'h',h,'rmsecv',0,'rmsep',0,'plots',0,'classic',classic,'rmsep_value',pred.rmsep,'rmsecv_value',rmsecv_value,'rsquared_value',R2,'rss_value',rss); 
    end
else
    if rmsecv
        error(['Both RMSECV and k were given.', ...
            'Please rerun your analysis with one of these inputs.  (see help file)'])
    end
%First stage: Obtain the scores T by first performing ROBPCA on z: 
    if st<=1
        out.robpca=robpca(z,'k',kr,'h',h,'plots',plotsrobpca,'kmax',kmaxr,'classic',classic,'mcd',mcd);
        out.h=h;
        out.centerz=out.robpca.M;
        out.sigmaxy=out.robpca.P(1:p1,:)*diag(out.robpca.L)*out.robpca.P(p1+1:p1+q1,:)';
        out.sigmax=out.robpca.P(1:p1,:)*diag(out.robpca.L)*out.robpca.P(1:p1,:)';
        out.xcentr=x-repmat(out.centerz(1:p1),n,1);
        out.ycentr=y-repmat(out.centerz(p1+1:p1+q1),n,1);
        out.weights2=out.robpca.flag.all; 
    end
    if st
        i=k;
    else 
        i=1;
    end
    while i<=k
        out.sigmayx=out.sigmaxy';
        if q1>p1        
            [RR,LL]=eig(out.sigmaxy*out.sigmayx); 
            [LL,I]=greatsort(diag(LL));
            rr=RR(:,I(1));
            qq=out.sigmayx*rr;  						                      
            qq=qq/norm(qq); 
        else
            [QQ,LL]=eig(out.sigmayx*out.sigmaxy);	
            [LL,I]=greatsort(diag(LL));
            qq=QQ(:,I(1));
            rr=out.sigmaxy*qq;
            rr=rr/norm(rr); 
        end
        tt=out.xcentr*rr;
        uu=out.ycentr*qq;	
        pp=out.sigmax*rr/(rr'*out.sigmax*rr);
        vv=pp;
        if i>1 												
            vv=vv-out.v*(out.v'*pp);
        end
        if vv'*vv==0
            error('The number of components is too large')
        end
        vv=vv./norm(vv);									
        out.sigmaxy=out.sigmaxy-vv*(vv'*out.sigmaxy); 
        out.v(:,i)=vv;
        out.q(:,i)=qq;
        out.t(:,i)=tt;
        out.u(:,i)=uu;
        out.p(:,i)=pp;
        out.r(:,i)=rr;
        i=i+1;
    end
    
    %Second Stage : Robust ROBPCA-regression
    robpcareg=robpcaregres(out.t,y,out.weights2);
    breg=robpcareg.coeffs(1:k,:);
    Yhat=out.t*breg+repmat(robpcareg.coeffs(k+1,:),n,1);
    b=out.r*breg;  
    int=robpcareg.coeffs(k+1,:)-out.centerz(1:p1)*out.r*breg;

    if rmsep==1
        rmse=rrmse(x,y,h,kmax,'RSIMPLS',0,k);  
        out.rmsep=rmse.rmsep;
    end
    
    % testing several output parameters
    if ~isnan(rmsep_value)
        out.rmsep = rmsep_value;
        options.rmsep = 1;
    end
    
    if ~isnan(rsquared_value)
        out.rcs = rsquared_value;
    end
    if ~isnan(rss_value)
        out.rcs = [out.rcs;sqrt(rss_value)];
    end
    if ~isnan(rmsecv_value)
        gammahalf = 0.5*sqrt(rss_value) + 0.5*rmsecv_value;
        out.rcs = [out.rcs;gammahalf;rmsecv_value];
        options.rmsecv = 1;
    end
    if any(isnan(rsquared_value)) && any(isnan(rss_value)) && any(isnan(rmsecv_value))
        out.rcs = 0;
    end

   %The output:
   out.T=out.t;
   out.weights.p=out.p;
   out.weights.r=out.r;
   out.kr=kr;
   out.h=h;
   out.alpha=alfa;
   out.slope=b;
   out.int=int;
   out.yhat=x*b+repmat(int,n,1);
   out.x=x;
   out.y=y;
   out.res=y-out.yhat;
   out.class='RSIMPLS';
   out.k=k;
   out.cov=robpcareg.cov;
    
   if ~st
       %calculation of robust distances
       %Score distance
       out.Tcov=robpcareg.sigma(1:k,1:k);
       out.Tcenter=robpcareg.center(1:k);
       out.sd=sqrt(mahalanobis(out.t,out.Tcenter,'cov',out.Tcov))';
       out.cutoff.sd=sqrt(chi2inv(0.975,k));
       %robust residual distance
       if q1==1
           out.resd=out.res/sqrt(out.cov);
       else
           out.resd=sqrt(mahalanobis(out.res,zeros(1,q1),'cov',out.cov))';
       end
       %robust orthogonal distances
       xtilde=out.t*out.p';
       Cdiff=out.xcentr-xtilde;
       for i=1:n
           out.od(i,1)=norm(Cdiff(i,:));
       end
       r=rank(x);
       if k~=r
           [m,s]=unimcd(out.od.^(2/3),out.h);
           out.cutoff.od = sqrt(norminv(0.975,m,s)^3);
       else
           out.cutoff.od=0;
       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);


       %Multivariate Rsquared
       Yw=y(out.flag.all==1,:);
       cYw=mcenter(Yw);
       res=out.res(out.flag.all==1,:);
       out.rsquared=1-(det(res'*res)/det(cYw'*cYw));

       %Assigning output
       if options.rmsep~=1 && options.rmsecv~=1
           out.rmsep=0;
       end

       if classic
           resultclassic=csimpls(x,y,'k',k,'plots',0);
       else
           resultclassic=0;
       end
       result=struct('slope',{out.slope}, 'int',{out.int},'fitted',{out.yhat},'res',{out.res}, 'cov',{out.cov},...
           'T',{out.T}, 'weights', {out.weights},'Tcenter',{out.Tcenter},'Tcov',{out.Tcov},'rsquared',{out.rsquared},'rcs',{out.rcs},'rmsep',{out.rmsep},...
           'k',{out.k},'alpha',{out.alpha},'h',{out.h},'sd', {out.sd},'od',{out.od},...
           'resd',{out.resd},'cutoff',{out.cutoff},'flag',{out.flag},'class',{out.class},'classic',{resultclassic},'robpca',{out.robpca});
       if result.rcs==0
           result=rmfield(result,'rcs');
       end
       if result.rmsep==0
           result=rmfield(result,'rmsep');
       end
   else
       result=out;
   end
end

% Plots
try
    if plots && options.classic
        makeplot(result,'classic',1,'labsd',labsd,'labod',labod,'labresd',labresd)
    elseif plots
        makeplot(result,'labsd',labsd,'labod',labod,'labresd',labresd)
    end
catch %output must be given even if plots are interrupted
end