clear all
cd ('\\stat.ad.ncsu.edu\Redirect\cwang19\Desktop\P3simulations')
addpath(genpath('./fslicing'))
addpath(genpath('./release2.16'))
rng(2016)
ntrain = 100;

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;

regular = 0; % regular =  0 means creating sparse data
nsim = 20;



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

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



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

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



& 0.7731 & 0.7769 &0.7836 &0.7895 &0.9507 &0.8950 &0.8354 &0.8141 &
& 0.0050 & 0.0052 &0.0051 &0.0052 &0.0049 &0.0053 &0.0047 &0.0049 &


& 0.8361 & 0.9586 &1.0122 &1.0149 &0.9836 &1.0359 &1.0184 &1.0184 &
& 0.0148 & 0.0094 &0.0014 &0.0000 &0.0085 &0.0015 &0.0003 &0.0003 &


& 0.8334 & 0.8261 &0.9249 &0.9644 &0.9954 &1.0357 &1.0167 &1.0167 &
& 0.0293 & 0.0417 &0.0310 &0.0225 &0.0175 &0.0031 &0.0004 &0.0004 &