function result = cvRobpca(data,kmax,resrob,rawres,h,csteps)

%CVROBPCA calculates the robust cross-validated PRESS (predicted residual error sum of squares) curve
% for ROBPCA in a fast way. This curve can be used to make a selection of the optimal number of 
% components. The function is used in robpca.m. 
%
% Input arguments:
%    data : the whole data set
%    kmax : the maximal number of components to be considered.
%  resrob : result of robpca for k = kmax on the whole data set. 
%  rawres : Optional input parameter. If equal to 1 (default), then the R-PRESS is calculated
%           based on the orthogonal distances of the points, that are calculated without performing the 
%           MCD step within ROBPCA. 
%       h : the quantile used in ROBPCA.
%  csteps : optional : csteps.value : whether c-steps should be performed (default) or not.                        
%                      csteps.number: the number of c-steps that must be performed. 
%
% Output: 
%   result.press    : the R-PRESS values when the minimum is used to define the weights.
%   result.weights  : the minimum, median and smooth weights.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at: 
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Sanne Engelen 
% Last Update: 01/07/2004, 03/07/2006
% Last Revision: 03/07/2006


n = size(data,1);
p = size(data,2);
r = rank(data);
alfa = 0.75;
teller_if_lus = 0;

% some initialisations
Pk_min_i = [];
Lk_min_i = [];
muk_min_i = [];
Tik = [];

if nargin < 4
    rawres = 1;
    alfa=0.75;
    h = min(floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*alfa),n);
    csteps.value = 1;
    csteps.number = 2;
elseif nargin == 4
    alfa=0.75;
    h = min(floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*alfa),n);
    csteps.value = 1;
    csteps.number = 2;
elseif nargin == 5
    csteps.value = 1;
    csteps.number = 2;
else
    csteps.value = csteps.value;
    if csteps.value
        csteps.number = csteps.number;
    end
end

outWeights = weightscvRobpca(data,r,resrob,kmax,h);
result.weights = outWeights;
w_min = outWeights.w_min;

% inputH0
H0 = resrob.Hsubsets.H0;
same.value = 0;

for i = 1:n
    same.value = 0;
    if isempty(find(H0 == i))
        if teller_if_lus >= 1
            same.value = 1;
        end
        teller_if_lus = teller_if_lus + 1;
    end
    
    if ~rawres
        % inputH1
        H1 = resrob.Hsubsets.H1;
        
        % Hfreq 
        Hfreq = resrob.Hsubsets.Hfreq;
    
        Hsets = [H0;H1;Hfreq];
        factor = 1;
        Hsets_min_i = RemoveObsHsets(Hsets,i);
        res = removeObsRobpca(data,i,kmax,Hsets_min_i,same,factor,csteps);
    else
        data_min_i = removal(data,i,0);
        H0_min_i = RemoveObsHsets(H0,i);
        if ~same.value
            [Pk_min_i,TH0,LH0_min_i,rH0_min_i,centerX,muk_min_i]=kernelEVD(data_min_i(H0_min_i,:));
            res.Pk_min_i = Pk_min_i;
            res.muk_min_i = muk_min_i;
        else
            res = same.res;
            Pk_min_i = res.Pk_min_i;
            muk_min_i = res.muk_min_i;
        end
    end

    if isempty(find(H0 == i))
        same.res = res;
    end
    
    Pkmax_min_i = res.Pk_min_i;
    muk_min_i = res.muk_min_i;
    
    for k = 1:kmax
        xhoed_k(i,(k-1)*p + 1:k*p) = (data(i,:) -  muk_min_i)*Pkmax_min_i(:,1:k)*Pkmax_min_i(:,1:k)' + muk_min_i;
        if k~=r
            odk(i,k) = norm(data(i,:) - xhoed_k(i,(k-1)*p + 1:k*p));
        else
            odk(i,k) = 0;
        end
    end
end

for k = 1:kmax
    press_min(k) = 1/sum(w_min)*w_min*odk(:,k).^2;
end

result.press = press_min;
%-----------------------------------------------------------------------------------------------------------
function Hsets_min_i = RemoveObsHsets(Hsets,i)

% removes the right index from the $h$-subsets in Hsets to 
% obtain (h - 1)-subsets.
% every h-subset is put as a row in Hsets.
% i is the index of the observation that is removed from the whole data.

for r = 1:size(Hsets,1)
    if ~isempty(find(Hsets(r,:)== i))
        Hsets_min_i(r,:) = removal(Hsets(r,:),0,find(Hsets(r,:) == i));
    else
        Hsets_min_i(r,:) = Hsets(r,1:(end-1));
    end

    for j = 1:length(Hsets_min_i(r,:))
        if Hsets_min_i(r,j) > i
            Hsets_min_i(r,j) = Hsets_min_i(r,j) - 1;
        end
    end
end
%-----------------------------------------------------------------------------------------------------------
function out = weightscvRobpca(data,r,resrob,kmax,h)

% computes the weights needed for the calculation of the R-PRESS.
%
% input:
%   data : the whole data set.
%      r : the rank of the data set.
%   resrob: the result of robpca on the whole data set for k = kmax.
%   kmax : the maximal number of components to be considered.
%   h : (n-h+1) measures the number of outliers the algorithm should 
%       resist. Any value between n/2 and n may be specified. 
%
% output: 
%   out.w_min    : the minimum weights


odk = [];
n = size(data,1);
p = size(data,2);

% Determine fixed weights:
for k = 1:kmax
    muk = resrob.M;
    xhoed_k(:,(k-1)*p + 1:k*p) = (data -  repmat(muk,n,1))*resrob.P(:,1:k)*resrob.P(:,1:k)' + repmat(muk,n,1);
    
    for i = 1:n
        % defining the od for each k
        if k~=r
            odk(i,k) = norm(data(i,:) - xhoed_k(i,(k-1)*p + 1:k*p));
        else
            odk(i,k) = 0;
        end
    end
    
    % defining weights for odk:
    if k~=r
        [m,s]=unimcd(odk(:,k).^(2/3),h);
        cutoff(k)=sqrt(norminv(0.975,m,s).^3);
        wod(:,k) = (odk(:,k) <= cutoff(k));
    else
        cutoff(k) = 0;
        wod(:,k) = 1;
    end
end

% determine the weights for every observation:
wk = wod;

if size(wk,1) == 1 || size(wk,2) == 1
    w_min = wk';
else
    w_min = min(wk,[],2)';
end

out.w_min = w_min;
out.odk = odk;
out.cutoff = cutoff;