clear;
addpath('../aux'); %test;
addpath('../figures');
addpath('../results'); %test;

design='NP'; %NP, SSG, HLLT
[f,sim,x_vis,f_vis]=get_design(design);

% sample size
N=1000;

% simulate data, split into stage 1 and 2 samples
rng('default');
[x,y,z]=sim(f,N);


%% kernel ridge reg - no IV, 2-fold CV

% precalculate all kernel matrices
df=get_K(x,y,z,x_vis);

% initialize hyperparameters for tuning
lambda0=log(0.05);

% tuning by CV
kernel_ridge_obj=@(lambda) kernel_ridge_cv(df,exp(lambda)); %exp to ensure pos
lambda_star=fminunc(kernel_ridge_obj,lambda0);

% evaluate on full sample using tuned hyperparameters
y_vis=kernel_ridge_pred(df,exp(lambda_star),3); %exp to ensure pos

% visualize estimator
plot(x_vis,f_vis,'LineWidth',5);
hold on;
scatter(x,y,36,[.5 .5 .5]);
plot(x_vis,y_vis,'--r','LineWidth',5);
xlabel('x','FontSize',20)
ylabel('y','FontSize',20)
legend({'Structural function','Data','KernelReg'},'Location','southeast','FontSize',20);
hold off;
hold off;
saveas(gcf,fullfile('../figures',strcat('kernel_ridge_',design)),'epsc');

%% kernel_ridge - 100 simulations

clear;
rng('default')
design='NP'; %NP, SSG, HLLT
[f,sim,x_vis,f_vis]=get_design(design);

% sample size
N=1000;
n_trials=100;
results=zeros(n_trials,N);
results_mse=zeros(n_trials,1);

for i=1:n_trials
    if mod(i,10)==0
        disp(num2str(i));
    end
    results(i,:)= sim_pred(design,N)';
    results_mse(i)=mse(results(i,:)',f_vis);
end

% visualize 100 samples
means=mean(results);
sorted=sort(results);
fifth=sorted(5,:);
nintetyfifth=sorted(95,:);

plot(x_vis,f_vis,'LineWidth',2);
hold on;
plot(x_vis,means,'--r','LineWidth',2);
plot(x_vis,fifth,':r','LineWidth',2);
plot(x_vis,nintetyfifth,':r','LineWidth',2);
xlabel('x');
ylabel('y');
hold off;
saveas(gcf,fullfile('../figures',strcat('kernel_ridge_',design,'_',num2str(n_trials))),'epsc');