function result=rpcr(x,y,varargin)

%RPCR is a 'Robust Principal Components Regression' method based on ROBPCA.
% 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. First, a robust principal components method is applied to the 
% predictor variables x (see robpca.m). Then a robust regression is performed. 
% For univariate y, the LTS regression is used (see ltsregres.m). When there are 
% several response variables, the MCD-regression method is applied (see mcdregres.m).
%
% The RPCR method is described in: 
%    Hubert, M., Verboven, S. (2003),
%    "A robust PCR method for high-dimensional regressors",
%    Journal of Chemometrics, 17, 438-452.
%
% 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 principal components to be used in the PCA step 
%                (default = min(rank(x),kmax)). If k is not specified,
%                it can be selected using the option 'rmsecv'. 
%         kmax : Maximal number of principal components to be used (default = 10).
%                If k is provided, kmax does not need to be specified, unless k is larger
%                than 10.                 
%        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 robust cross-validated 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.
%    intadjust : If equal to one, the intercept adjustment for the
%                LTS-regression will be calculated (default = 0). See ltsregres.m for
%                details.
%        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 the input argument 'classic' is equal to one, the classical
%                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 diagnostic plots. (default = 3)
%        labod : The 'labod' observations with largest orthogonal distance are
%                labeled on the score outlier map. (default = 3)    
%      labresd : The 'labresd' observations with largest residual distance are
%                labeled on the regression outlier map. (default = 3)   
%      classic : If equal to one, the classical PCR analysis will be performed as well
%                (see also cpcr.m). (default = 0)
% 
%
% I/O: result=rpcr(x,y,'k',0,'kmax',10,'alpha',0.75,'h',h,'rmsecv',0,'rmsep',0,...
%             'plots',1,'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: result=rpcr(x,y,'alpha',0.65,'k',5,'plots',0,'classic',1);
%          result=rpcr(x,y,'kmax',5,'labresd',7);
%
% The output of RPCR is a structure containing:
%
%   result.slope     : Robust slope
%   result.int       : Robust intercept 
%   result.fitted    : Robust prediction vector
%   result.res       : Robust residuals
%   result.cov       : Estimated variance-covariance matrix of the residuals  
%   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 sfor
%                      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 principal components
%   result.h         : The quantile h used throughout the algorithm
%   result.sd        : Robust score distances within the robust PCA subspace
%   result.od        : Robust orthogonal distances to the robust PCA subspace 
%   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 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     : 'RPCR'
%   result.classic   : If the input argument 'classic' is equal to one, this structure
%                      contains results of the classical PCR analysis (see also cpcr.m). 
%   result.robpca    : Full output of the robust PCA analysis (see robpca.m) 
%   result.lts       : If there is one response variable: full output of the LTS regression.
%                      (see ltsregres.m)
%   result.mcdreg    : If there are several response variables: full output of the MCD regression.
%                      (see mcdregres.m)
%
% 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
% Created on: 31/10/2000
% Last Update: 09/04/2004,  03/07/2006
% Last Revision: 04/08/2006

%
% checking input
%
if rem(nargin,2)~=0
    error('Number of inputarguments must be even!');
end
%
% initialization with defaults
%
counter=1;
[n,p]=size(x);
r=rank(x);
[n2,q]=size(y);
if n~=n2
    error('The response variables and the predictor variables have a different number of observations.')
end
kmax=min([10,floor(n/2),r]);
alfa=0.75;
k=0;
h=floor(2*floor((n+kmax+2)/2)-n+2*(n-floor((n+kmax+2)/2))*alfa);
default=struct('intadjust', 0, 'alpha',alfa,'h',h,'k',k,'kmax',kmax,'plots',1,...
    'rmsecv',0,'classic',0,'labsd',3,'labod',3,'labresd',3,'rmsep',0);
list=fieldnames(default);
options=default;
IN=length(list);
i=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 
    dummy=sum(strcmp(chklist,'h')+2*strcmp(chklist,'alpha'));
    switch dummy
        case 0 % defaultvalues should be taken
            options.alpha=alfa;
            options.h=h;
        case 3
            error('Both inputarguments 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');
        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
    options.h=floor(options.h);
    options.kmax=floor(options.kmax);
    options.k=floor(options.k);
    options.labsd=max(0,min(floor(options.labsd),n));
    options.labod=max(0,min(floor(options.labod),n));
    options.labresd=max(0,min(floor(options.labresd),n));
    kmax=min([options.kmax,floor(n/2),r]);
    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
    if options.k < 0 
        options.k=0;
    elseif options.k > kmax 
        mess=sprintf(['Attention (rpcr.m): The required number of principal components, k = ',num2str(options.k)...
            ,'\n is larger than kmax= ',num2str(kmax),'; k is set to ',num2str(kmax)]);
        disp(mess)
        options.k=kmax;
    end
    if dummy==1% checking inputvariable h
        if options.h-floor(options.h)~=0
            mess=sprintf(['Attention (rpcr.m): h must be an integer. \n']);
            disp(mess)
        end
        if options.k==0
            if options.h<floor((n+kmax+1)/2 )
                options.h = quanf(0.5,n,kmax);
                mess=sprintf(['Attention (rpcr.m): h should be larger than (n+2)/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+options.k+1)/2) % checking inputvariable h
                options.h = quanf(0.5,n,k);
                mess=sprintf(['Attention (rpcr.m): h should be larger than (n+2)/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 options.k==0
                options.h=floor(2*floor((n+kmax+2)/2)-n+2*(n-floor((n+kmax+2)/2))*options.alpha);
            else
                options.h=floor(2*floor((n+options.k+1)/2)-n+2*(n-floor((n+options.k+1)/2))*options.alpha);%k+1 ipv k?
            end
            mess=sprintf(['Attention (rpcr.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 (rpcr.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 (rpcr.m): Alpha should be smaller than 1.\n',...
                    'It is set to 0.75.']);
            disp(mess)
        end
        if options.k==0
            options.h=floor(2*floor((n+kmax+2)/2)-n+2*(n-floor((n+kmax+2)/2))*options.alpha);
        else
            options.h=floor(2*floor((n+options.k+1)/2)-n+2*(n-floor((n+options.k+1)/2))*options.alpha);%k+1 ipv k?
        end
    end
end
%
%
% MAIN PART
%
% If y is univariate, perform LTS regression, else MCD regression
if q==1 
    options.method='lts';
else
    options.method='mcd';
end
% Default value for the number of latent variables k
if n<kmax+2
     error('You need more data points to compute this estimator.')
end
%  
% initializing and checking k
if options.k==0 % no optimal number of latent variables given by user
    kfixed = 0;
    if options.rmsecv==0 % no cross validation asked -> calculate R2
        [R2,final]=rsquared(x,y,kmax,'RPCR',options.h);
        rss = final.rss;
        options.k=final.k;
        outrobpca=robpca(x,'h',options.h,'k',options.k,'kmax',kmax,'plots',0,'classic',options.classic);
        options.k=outrobpca.k;
        final=regression(x,y,options.k,outrobpca,options);
        if options.rmsep==1 % the rmsep is asked, without prior knowledge
            rmse=rrmse(x,y,options.h,kmax,'RPCR',0,options.k);
            final.rmsep=rmse.rmsep;
        end
        final.robpca=outrobpca;
        final.R2 = R2;
        final.rss = rss;
    else % RMSE with cross validation
        crosscv=rrmse(x,y,options.h,kmax,'RPCR'); 
        options.k=crosscv.k;
        rmsecv=crosscv.rmsecv;
        pred=rrmse(x,y,options.h,kmax,'RPCR',0,options.k,crosscv.weight,crosscv.res); %calculate rmsep with rmsecv weights and residuals
        outrobpca=robpca(x,'h',options.h,'k',options.k,'kmax',kmax,'plots',0,'classic',options.classic);
        options.k=outrobpca.k;
        final=regression(x,y,options.k,outrobpca,options);
        final.robpca=outrobpca;
        final.rmsecv=rmsecv;
        final.rmsep=pred.rmsep;
        final.R2 = crosscv.R2;
        final.rss = crosscv.rss;
    end
else % optimal number of latent variables is given by user
    if options.rmsecv
        error(['Both RMSECV and k were given.', ...
            'Please rerun your analysis with one of these inputs.  (see help file)'])
    end
    kfixed = 1;
    k=min(r,min(options.k,kmax));
    options.k=k;
    outrobpca=robpca(x,'h',options.h,'k',options.k,'kmax',kmax,'plots',0,'classic',options.classic);
    options.k=outrobpca.k;
    final=regression(x,y,options.k,outrobpca,options);
    if options.rmsep==1 % the rmsep is asked, without prior knowledge
        rmse=rrmse(x,y,options.h,kmax,'RPCR',0,options.k);
        final.rmsep=rmse.rmsep;
    end
    final.robpca=outrobpca;
end

final.h=options.h;
final.k=options.k;
final.alpha=options.alpha;

% scores-distances
final.sd=final.robpca.sd;
quan=chi2inv(0.975,final.k);
final.cutoff.sd=quan^0.5;

% orthogonal distances 
final.od=final.robpca.od;
final.cutoff.od=final.robpca.cutoff.od;
final.x=x;
final.y=y;

if options.classic==1
    final.classic=cpcr(x,y,'k',final.k,'plots',0);
else
    final.classic=0;
end

if strmatch(options.method,'mcd','exact') 
    final.resd=final.mcdreg.resd;
    final.cutoff.resd=final.mcdreg.cutoff.resd;
else
    final.resd=final.lts.res/final.lts.scale;
    final.cutoff.resd=sqrt(chi2inv(0.975,1));
end
final.class='RPCR';

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

% Assigning output
if kfixed == 0
    final.rcs = [final.R2;sqrt(final.rss)];
else
    final.rcs = 0;
end

if options.rmsecv~=0
    gammahalf = 0.5*sqrt(final.rss) + 0.5*final.rmsecv;
    final.rcs = [final.rcs;gammahalf;final.rmsecv];
    rmsecv_ind = 1;
else
    rmsecv_ind = 0;
end
if options.rmsep~=1 && rmsecv_ind == 0
    final.rmsep=0;
end

if isfield(final,'lts')
    result=struct('slope',{final.slope}, 'int',{final.int}, 'fitted',{final.fitted},'res',{final.res},...
        'cov',{final.cov},'rsquared',{final.rsquared},'rcs',{final.rcs},'rmsep',...
        {final.rmsep},'k',{final.k},'alpha',{final.alpha},'h',{final.h},'sd', {final.sd},'od',{final.od},'resd',{final.resd},...
        'cutoff',{final.cutoff},'flag',{final.flag},'class',{final.class},'classic',{final.classic},...
        'robpca',{final.robpca},'lts',{final.lts});
else
    result=struct('slope',{final.slope}, 'int',{final.int}, 'fitted',{final.fitted},'res',{final.res},...
        'cov',{final.cov},'rsquared',{final.rsquared},'rcs',{final.rcs},'rmsep',{final.rmsep},...
        'k',{final.k},'alpha',{final.alpha},'h',{final.h},'sd', {final.sd},'od',{final.od},'resd',{final.resd},...
        'cutoff',{final.cutoff},'flag',{final.flag},'class',{final.class},'classic',{final.classic},...
        'robpca',{final.robpca},'mcdreg',{final.mcdreg});
end
if result.rcs == 0
    result = rmfield(result,'rcs');
end
if result.rmsep==0 
    result=rmfield(result,'rmsep');
end

% Plots
try
    if options.plots && options.classic
        makeplot(result,'classic',1,'labsd',options.labsd,'labod',options.labod,'labresd',options.labresd)
    elseif options.plots
        makeplot(result,'labsd',options.labsd,'labod',options.labod,'labresd',options.labresd) 
    end
catch %output must be given even if plots are interrupted 
    %> delete(gcf) to get rid of the menu 
end
%---------------------------------
function out=regression(x,y,d,pre_out,options)

switch options.method
case 'mcd' % MCD regression for multivariate responses
    out.mcdreg=mcdregres(pre_out.T(:,1:d),y,'k',options.k,'h',options.h,'plots',0);
    if ~isstruct(out.mcdreg)
        return
    end
    %%coefficients in the original space
    out.slope=pre_out.P(:,1:d)*out.mcdreg.slope;
    out.int=out.mcdreg.int-pre_out.M*pre_out.P(:,1:d)*out.mcdreg.slope;
    out.fitted=[pre_out.T(:,1:d) ones(size(pre_out.T(:,1:d),1),1)]*[out.mcdreg.slope; out.mcdreg.int]; 
    out.res=y-out.fitted;
    out.cov=out.mcdreg.cov; % covariance matrix of residuals (capital sigma)
    out.name='Multivariate robust MCD-regression';  
    out.rsquared=out.mcdreg.rsquared;
case 'lts' 
    [out.lts,raw]=ltsregres(pre_out.T(:,1:d),y,'plots',0,'h',options.h,'intadjust',options.intadjust);
    %%coefficients in the original space;
    out.lts.weights=raw.wt;
    out.slope=pre_out.P(:,1:d)*out.lts.slope;
    out.int=out.lts.int-pre_out.M*pre_out.P(:,1:d)*out.lts.slope;
    out.fitted=[pre_out.T(:,1:d) ones(size(pre_out.T(:,1:d),1),1)]*[out.lts.slope; out.lts.int];
    out.res=y-out.fitted; 
    out.cov=out.lts.scale.^2;
    if out.cov==0
        mess=sprintf(['Attention (rpcr.m): No standardized residuals could be calculated \n',...
            'because their robust scale is zero.']);
        disp(mess)
        out.stdres=NaN;
    else
        out.stdres=out.res/out.lts.scale;
    end
    out.rsquared=out.lts.rsquared;
    out.name='Robust LTS regression';
end

%-----------------------------------------------------------------------------------------
function quan=quanf(alfa,n,rk)

quan=floor(2*floor((n+rk+1)/2)-n+2*(n-floor((n+rk+1)/2))*alfa);