clear all
cd ('\Users\chongwang\Documents\SAS backup\P3simulations')
cd \\stat.ad.ncsu.edu\Redirect\cwang19\Desktop\P3simulations
addpath(genpath('./fslicing'))
addpath(genpath('./release2.16'))
rng(2016)
ntrain = 400;

kn_all = [5 10 20 30];  %sn=kn:  how many eigenvalues kept in xcov


normtemp2 = 0;
normtemp1 = 0;

vary = 1;
varx = 0;

lintx = 10;
ngrid = 101;


regular = 2; % regular =  0 means creating sparse data
nsim = 100;
evals = sort((0.01:0.01:1),'descend'); C1 = 0.82; model = 1;  
uu_all = [1 2 3 4];    %uu: envelope dimension
dd_all = [2 2 2 2];    %dd: number of terms (powers) in the Krylov series


results = zeros(nsim,6);
for ss = 1:nsim
    % Create training sample
    [tx, x, y, t, ~, ~, truebeta] = getFunData_eg1(ntrain, lintx, ngrid, model, ...
regular, varx, vary, evals, C1);
    
    dimCS = size(truebeta,1);   %dimension of the central subspace
    [~, Ptruephi, ~, ~] = getEigens(truebeta'*truebeta,t,t,dimCS);
    Ptrue = Ptruephi*Ptruephi';
    
%     % Create testing sample
%     [t_x, newx, ~, ~, Ex, Ey, ~] = getFunData_eg1(ntest, lintx, ngrid, model, ...
%         regular, varx, vary, evals);
    
    for tuning=1:4  %three tuning parameter options for each of CUME and ECUME
        uu=uu_all(tuning); 
        dd=dd_all(tuning);
        kn = kn_all(tuning);
        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, betaJ, Jprojx, Jf_projx,...
            Jintercept,  betaPsi, PSIprojx, PSIf_projx,...
            PSIintercept, edrout, Jgamout,...
            PSIgamout]= edrIndex(x, tx, y, p);
        
        %%%%%%%%%%%%%%%%%%%%%%
        %CUME                %
        %%%%%%%%%%%%%%%%%%%%%%
%         [~, Phatphi, ~, ~] = getEigens(betaJ*betaJ', t ,t_beta, dimCS);
%         Phat = Phatphi*Phatphi';
%         normtemp1 = mean(svds(Ptrue-Phat, 2*dimCS)*range(t)/(length(t)-1));
%         clear next newprojx newf_projx newy Phatphi Phat ind
        
        %%%%%%%%%%%%%%%%%%%%%%
        %Envelop CUME        %
        %%%%%%%%%%%%%%%%%%%%%%
        [~, Phatphi, ~, ~] = getEigens(betaPsi*betaPsi', t ,t_beta, dimCS);
        Phat = Phatphi*Phatphi';
        normtemp2 = mean(svds(Ptrue-Phat, 2*dimCS)*range(t)/(length(t)-1));
        clear next newprojx newf_projx newy Phatphi Phat ind
        
        results(ss,[tuning 4+tuning]) = [normtemp2 normtemp1];
    end
    
end

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








