function [cvuu, cvdim] = cv2(n, folds, alltx, allx, ally, t, allEx, allEy, ngrid)

%input the whold traing data set
%do a ten folds cross validation
%output: selected dimension 
kn_all = [5 10 20 30];  %sn=kn:  how many eigenvalues kept in xcov
uu_all = [1 2 3 4];    %uu: envelope dimension   
regular = 2;
Indices = crossvalind('Kfold', n, folds);
results = zeros(folds ,10);
for indx = 1:folds
    j=0;
    
    % Create testing sample
    t_x = alltx(Indices == indx);
    Ex = allEx(Indices == indx);
    Ey = allEy(Indices == indx);
    yvar =  var(Ey);
    % Create training sample
    tx = alltx(Indices ~= indx);
    x = allx(Indices ~= indx);
    y = ally(Indices ~= indx);
    
    
    for tuning=1:4  %three tuning parameter options for each of CUME and ECUME
        uu=uu_all(tuning); 
        dd=uu;
        kn = kn_all(tuning);
        sdr_all = [1:uu];
        for i = 1:length(sdr_all)
            dimCS=sdr_all(i); 
            
            p = edrIndexOptions('regular', regular,   'verbose', 'off', 'ngrid', ngrid, 'kernel', 'gauss', ...
            'newdata', t, 'sn', kn, 'K', dimCS, 'edrMethod', 'Ecume', 'uu', uu , 'dd', dd , 'npoly', 1);
        
            [t_beta, ~, ~,betaPsi, PSIprojx, ~] = NEWedrIndex(x, tx, y, p);

            %%%%%%%%%%%%%%%%%%%%%%
            %Envelop CUME        %
            %%%%%%%%%%%%%%%%%%%%%%
            % Prediction
            %[next, newprojx, newf_projx, newy] = edrIndexPred('gaussian', t_beta, betaPsi, PSIintercept, PSIprojx, PSIf_projx, newx, t_x); 
            %[next, newprojx, newy] = NEW2edrIndexPred(t_beta, betaPsi, PSIprojx, y, Ex, t, criteria); 
            [next, ~, newy] = NEW2edrIndexPred(t_beta, betaPsi, PSIprojx, y, Ex, t_x); 
            ind = sum(next,1) == size(next,1); % Only consider non-extrapolated values ????????????????????????????
            % Calculate relative prediction error
            petemp2 = sum((Ey(ind) - newy(ind)).^2)/(sum(ind)*yvar);
            clear next newprojx newf_projx newy Phatphi Phat ind
            j=j+1;
            results(indx,j) = [petemp2];  
        end
        
    end
    
end

combination = [1 2 2 3 3 3 4 4 4 4; 1 1 2 1 2 3 1 2 3 4];
[~, ind] = min(mean(results),[],2);
cvuu = combination(1,ind);
cvdim = combination(2,ind);




