function res = removeObsRobpca(data,i,k,Hsets_min_i,same,factor_ind,csteps)

%REMOVEOBSROBPCA is an auxiliary function to perform cross-validation with ROBPCA, 
% RPCR, RSIMPLS, RSIMCA (see cvRobpca.m, cvRpcr.m, cvRsimpls.m).
%
% Input: 
%    data : the data set
%       i : the observation that is removed, index with respect to the whole data set.
%       k : the number of principal components that has to be calculated.
%  Hsets_min_i : contains H0_min_i, H1_min_i and Hfreq_min_i as first, second and third row respectively.
%                The h-subsets are implemented by means of indices of the observations in data_min_i, which is
%                is the original data set minus sample i. 
%  same : structure : 
%           same.value : indicates whether some part of this algorithm can be skipped (= 1) or not ( = 0). 
%           same.res : if same.value = 1, then some additional information is needed. 
%  factor : optional. default  = 0, then the original consistency factor is used.
%             Else factor = 1, the consistency factor is adapted to the kmax approach.
%  csteps : optional structure: 
%            csteps.value : 1 (default, then csteps are performed within robpca), 0 (no csteps are performed within robpca)
%            csteps.number : the number of csteps that need to be performed (default = 2).
%
% Output:
% res is the result structure. It contains:
%   Pk_min_i  : update of the loadingmatrix for a certain k when observation i is deleted.
%   Lk_min_i  : update of the eigenvalues for a certain k when observation i is deleted.
%   muk_min_i : update of the center for a certain k when observation i is deleted.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by S. Engelen

rot = [];
center = [];
P1 = [];

n = size(data,1);
p = size(data,2);
h = size(Hsets_min_i,2) + 1;

if nargin < 6    
    factor_ind = 0;    
    csteps.number = 2;   
    csteps.value = 1;
elseif nargin == 6    
    csteps.number = 2;    
    csteps.value = 1;
end

data_min_i= removal(data,i,0);
mu0 = mean(data);
center = mu0;

H0_min_j = Hsets_min_i(1,:);
H1_min_t = Hsets_min_i(2,:);

if ~same.value   
    [PH0_min_i,TH0,LH0_min_i,rH0_min_i,centerX,mu1trafo]=kernelEVD(data_min_i(H0_min_j,:));  
    LH0_min_i = diag(LH0_min_i);   
    
    % calculation of the projection
    % res.T2tilde = T2tilde;   
    center = mu1trafo;    
    rot = PH0_min_i;    
    T2tilde = (data_min_i - repmat(mu1trafo,n-1,1))*PH0_min_i;  
    T2tilde = T2tilde(:,1:k);  
    rot = rot(:,1:k);       
    
    % defining the outputstructure res:  
    res.PH0_min_i = PH0_min_i; 
    res.LH0_min_i = LH0_min_i;   
    res.mu1trafo = mu1trafo;
else
    % defining the input structure input
    res = same.res;  
    PH0_min_i = res.PH0_min_i;  
    LH0_min_i = res.LH0_min_i;   
    mu1trafo = res.mu1trafo;    
    T2tilde = (data_min_i - repmat(mu1trafo,n-1,1))*PH0_min_i;  
    center = mu1trafo;    
    rot = PH0_min_i;  
    rot = rot(:,1:k);   
    T2tilde = T2tilde(:,1:k);
end

mah = libra_mahalanobis(T2tilde,zeros(1,k),'invcov',1./diag(LH0_min_i(1:k,1:k)));
oldobj = prod(diag(LH0_min_i(1:k,1:k)));P4 = eye(k);

if csteps.value    
    oldobj = prod(diag(LH0_min_i(1:k,1:k)));  
    for j = 1:csteps.number     
        [mahsort, indsort] = sort(mah);      
        dataH1 = T2tilde(indsort(1:(h-1)),:);    
        [P,T,L,r3,Xm,clmX] = classSVD(dataH1);  
        obj = prod(L);      
        T2tilde = (T2tilde - repmat(clmX,n-1,1))*P;    
        center = center + clmX*rot';   
        rot = rot*P;      
        mah = libra_mahalanobis(T2tilde,zeros(1,size(T2tilde,2)),'invcov',1./L);      
        P4 = P4*P;    
        if abs(oldobj - obj) > 1.e-12      
            oldobj = obj;    
        else
            break   
        end
    end
else
    obj = oldobj;
end

% extra reweighting:
if k~=r3  
    XRc= data_min_i-repmat(center,n-1,1);   
    Xtilde = T2tilde*rot'; 
    Rdiff = XRc-Xtilde;   
    for i=1:n     
        odh(i,1)=norm(Rdiff(i,:));   
    end
    [m,s]=unimcd(odh.^(2/3),h);    
    cutoffodh = sqrt(norminv(0.975,m,s).^3);     
    indexset = find(odh<=cutoffodh)';
    [P,Threw,Lrew,rrew,Xmrew,clmX]=kernelEVD(data_min_i(indexset,:));   
    center = clmX;   
    rot = P(:,1:k);
end
T2tilde = (data_min_i - repmat(center,n-1,1))*rot;

% Perform mcdcov on some H-subsets:
[res_min_i,raw_min_i] = mcdcov(T2tilde,'Hsets',Hsets_min_i,'h',h-1,'ntrial',250,'plots',0,'factor',1);
% perform the last part of ROBPCA:
if raw_min_i.objective < obj  
    z = res_min_i;
else
    sortmah = sort(mah);  
    if h==n      
        factor=1;  
    else
        if factor_ind == 1        
            factor = sortmah(h-1)/chi2inv((h-1)/size(T2tilde,1),size(T2tilde,2)/2); % adjusted factor.       
        else
            factor = sortmah(h-1)/chi2inv((h-1)/size(T2tilde,1),size(T2tilde,2)); % non adjusted factor    
        end
    end
    mah = mah/factor;  
    weights = mah <= chi2inv(0.975,size(T2tilde,2)); 
    [center_noMCD,cov_noMCD] = weightmecov(T2tilde,weights);  
    mah = libra_mahalanobis(T2tilde,center_noMCD,'cov',cov_noMCD);  
    z.flag = (mah <= chi2inv(0.975,size(T2tilde,2)));   
    z.center = center_noMCD;  
    z.cov = cov_noMCD;
end

covf=z.cov;
centerf=z.center;

% The final PC:
[P6tilde,L6]=eig(covf);
[L6,I]=greatsort(diag(L6));
P6tilde=P6tilde(:,I);
T_min_i=(T2tilde-repmat(centerf,n-1,1))*P6tilde;

Pk_min_i=rot(:,1:k)*P6tilde;
Lk_min_i = L6;

res.Pk_min_i = Pk_min_i;
res.Lk_min_i = Lk_min_i;
res.muk_min_i = center + centerf*rot';