clear all
cd ('\\stat.ad.ncsu.edu\Redirect\cwang19\Desktop\P3simulations')
addpath(genpath('./fslicing'))
addpath(genpath('./release2.16'))
addpath(genpath('./mvkernelsmoothing'))


kn_all = [5 10 20 30];  %sn=kn:  how many eigenvalues kept in xcov
uu_all = [1 2 3 4];    %uu: envelope dimension    
%dd_all = [1 1 1 1];    %dd: number of terms (powers) in the Krylov series
%dd_all = [1 2 2 2]; 


vary = 1;
varx = 0;
lintx = 10;
ngrid = 101;

%evals = sort((0.01:0.01:1),'descend'); C1 = 0.82; model = 1;  %model 1 is I;  %%%%%%%%%%%%%%%%%%%%%%(a)
%evals = sort([2.1:0.01:2.15 1./(7:100).^1.25],'descend');C1 = 0.55; model = 1;  %model 1 is I;  %%%(b)
%evals = sort([1.55:0.05:2.0 10./(11:100)],'descend'); C1 = 0.58; model = 1;  %model 1 is I;   %%%%%%%%%%%(c)

%evals = sort((0.01:0.01:1),'descend'); C1 = 3; model = 2;  %model 1 is I;  %%%%%%%%%%%%%%%%%%%%%%(a)
%evals = sort([2.1:0.01:2.15 1./(7:100).^1.25],'descend');C1 = 2.4; model = 2;  %model 1 is I;  %%%(b)
%evals = sort([1.55:0.05:2.0 10./(11:100)],'descend'); C1 = 2.4; model = 2;  %model 1 is I;   %%%%%%%%%%%(c)

regular = 2; % regular =  0 means creating sparse data
rng(2016)
n = 400;
folds = 10;
nsim = 10;

% Create data set
[alltx, allx, ally, t, allEx, allEy, ~] = getFunData_eg1(n, lintx, ngrid, model, ...
regular, varx, vary, evals, C1);
    


Indices = crossvalind('Kfold', n, 10);


combination = [1 2 2 3 3 3 4 4 4 4; 1 1 2 1 2 3 1 2 3 4];

results = zeros(nsim ,10);
for ss = 1:nsim 
    j=0;
    
    % Create testing sample
    t_x = alltx(Indices == ss);
    Ex = allEx(Indices == ss);
    Ey = allEy(Indices == ss);
    yvar =  var(Ey);
    % Create training sample
    tx = alltx(Indices ~= ss);
    x = allx(Indices ~= ss);
    y = ally(Indices ~= ss);
    
    
    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', 101, 'kernel', 'gauss', ...
            'newdata', t, 'sn', kn, 'K', dimCS, 'edrMethod', 'Ecume', 'uu', uu , 'dd', dd , 'npoly', 1);
        
            [t_beta, ~, ~,betaPsi, PSIprojx, edrout] = 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, newprojx, newy] = NEWedrIndexPred(t_beta, betaPsi, PSIprojx, y, Ex, t); 
            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(ss,j) = [petemp2];  
        end
        
    end
    
end

%mean(results);
%std(results)./sqrt(nsim);

[val, ind] = min(results,[],2)
values = [1:10]
count = histc(ind,values)


formatSpec = '& %4.4f & %4.4f &%4.4f &%4.4f &%4.4f &%4.4f &%4.4f &%4.4f &%4.4f &%4.4f &\n';
formatSpec0 = '& %4.0f & %4.0f &%4.0f &%4.0f &%4.0f &%4.0f &%4.0f &%4.0f &%4.0f &%4.0f &\n';
fprintf(formatSpec0,count)
fprintf(formatSpec,mean(results))
fprintf(formatSpec,std(results)./sqrt(nsim ))

plot(mean(results))



