%% clear all variables and start fresh
clear variables

%% set seed
seed = 6;
rng(seed)

%% problem parameters
n = 2500;
p = 2000;

n_test = 1e5;

%% generate train and test features

feat_cov_type = 'iso';
feat_cov = generate_feature_covariance(p, feat_cov_type);

feat_type = 't';
feat_param = 5;

X = generate_feature_matrix(n, p, feat_cov, feat_type, feat_param);
X_test = generate_feature_matrix(n_test, p, feat_cov, feat_type, feat_param);

%% generate train and test responses

% signal generation
sig_type = 'iso';
sig_energy = 0.5;
beta0 = sqrt(sig_energy) * generate_signal_vector(p, sig_type);

% noise generation
noi_type = 't';
noi_param = 5;
noi_energy = 1;
eps = sqrt(noi_energy) * generate_noise_vector(n, noi_type, noi_param);
eps_test = sqrt(noi_energy) * generate_noise_vector(n_test, noi_type, noi_param);

% response generation
mod_type = 'nlsq';
y = generate_response_vector(X, beta0, eps, mod_type, feat_cov);
y_test = generate_response_vector(X_test, beta0, eps_test, mod_type, feat_cov);

%% fit estimator and get smoother

lambda = 0;
ridge_estim = fit_ridge(X, y, lambda);

ridge_smoother = get_ridge_smoother(X, lambda);

%% get predicted values and compute residuals

hy = predict_linestim(ridge_estim, X);
resid_orig = y - hy;

hy_test = predict_linestim(ridge_estim, X_test);
resid_test = y_test - hy_test;

resid_gcv = resid_orig / (1 - trace(ridge_smoother) / n);
resid_loo = resid_orig ./ (1 - diag(ridge_smoother));

if (lambda == 0) && (p > n)
   
    gram_mat = X * X' / n;
    resid_mod = gram_mat \ y;
    
    resid_gcv = resid_mod / (trace(pinv(gram_mat)) / n);
    resid_loo = resid_mod ./ diag(pinv(gram_mat));
    
end

%% plot histograms of residuals

h_oos = histogram(resid_test, ...
    'Normalization', 'pdf', ...
    'DisplayStyle', 'stairs', ...
    'LineWidth', 2, ...
    'DisplayName', 'Test');
hold on
grid on
% h_oos.BinWidth = h_gcv.BinWidth;
h_raw = histogram(resid_orig, ...
    'Normalization', 'pdf', ...
    'DisplayStyle', 'stairs', ...
    'LineWidth', 2, ...
    'DisplayName', 'Train');
% legend('oos', 'gcv', 'raw')
% h_gcv = histogram(resid_gcv, ...
%     'Normalization', 'pdf', ...
%     'DisplayStyle', 'stairs', ...
%     'LineWidth', 2, ...
%     'DisplayName', 'GCV');
h_gcv = histogram(resid_loo, ...
    'Normalization', 'pdf', ...
    'DisplayStyle', 'stairs', ...
    'EdgeColor', [0.4940 0.1840 0.5560], ...
    'LineWidth', 2, ...
    'DisplayName', 'LOOCV');
xlim([-10 10])
% ylim([0 1])
xlabel('Error')
ylabel('PDF')
legend show
legend boxoff
hold off

%% let's try nhist function

% histcell.oos = resid_test;
% histcell.gcv = resid_gcv;
% histcell.raw = resid_orig;
% 
% nhist(histcell, ...
%     'factor', 0.9, ...
%     'noerror')
% hold on
% grid on
% xlim([-4 4])
% ylim([0 1])

%
%% save the histogram plot

% if feat_param == 0.5
%     feat_param_str = 'h';
% elseif feat_param == 5
%     feat_param_str = num2str(feat_param);
% end
% 
% if noi_param == 0.5
%     noi_param_str = 'h';
% elseif noi_param == 0.8
%     noi_param_str = 'r';
% elseif noi_param == 0.2
%     noi_param_str = 'l';
% elseif noi_param == 5
%     noi_param_str = num2str(noi_param);
% end
% 
% if lambda < 1
%     lam_str = 'low';
% else
%     lam_str = num2str(lambda);
% end

save_str = strcat(...
    'hist_comp_loo', ...
    'n', num2str(n), '_', ...
    'p', num2str(p), '_', ...
    'nt', num2str(n_test), '_', ...
    'featcov', feat_cov_type, '_', ...
    'feat', feat_type, '_', ...
    'featparam', num2str(feat_param), '_', ...
    'sig', sig_type, '_', ...
    'noi', noi_type, '_', ...
    'noiparam', num2str(noi_param), '_', ...
    'mod', mod_type, '_', ...
    'lam', num2str(lambda), '_', ...
    'seed', num2str(seed));
% 
% save_folder = '../figures/';
% 
% saveas(h_oos, fullfile(save_folder, save_str), 'epsc')

saveas(gcf, save_str, 'epsc')
saveas(gcf, save_str, 'png')

%% plot cdfs of residuals

ecdf(resid_test)
hold on
grid on
ecdf(resid_orig)
ecdf(resid_loo)

legend('Test', 'Train', 'LOOCV')

xlabel('Error')
ylabel('PDF')
hold off

%% save the cdf comparison plot

save_str = strcat(...
    'cdf_comp_loo', ...
    'n', num2str(n), '_', ...
    'p', num2str(p), '_', ...
    'nt', num2str(n_test), '_', ...
    'featcov', feat_cov_type, '_', ...
    'feat', feat_type, '_', ...
    'featparam', num2str(feat_param), '_', ...
    'sig', sig_type, '_', ...
    'noi', noi_type, '_', ...
    'noiparam', num2str(noi_param), ...
    'mod', mod_type, '_', ...
    'lam', num2str(lambda), '_', ...
    'seed', num2str(seed));
% 
% save_folder = '../figures/';
% 
% saveas(h_oos, fullfile(save_folder, save_str), 'epsc')

saveas(gcf, save_str, 'epsc')
saveas(gcf, save_str, 'pdf')