rf_test = load('rainfall_for_risk.mat');
rf_stoc = load('ypred_stochastic.mat');
rf_corr = load('rainfall_for_risk_corrected.mat');

x0 = 0.001:0.001:0.999; x = 1./(1-x0);
year_range = 2010:2019;

y_pred = NaN(length(x0), length(year_range));
y_true = NaN(length(x0), length(year_range));
y_pred_stoc = NaN(length(x0), length(year_range));
y_pred_corr = NaN(length(x0), length(year_range));

for year=year_range
    p0 = rf_test.rf.(strcat('pred_rf_gp_gan_topo_', num2str(year)))(:);
    p1 = p0(p0>0);
    a = gpfit(p1);
    y_pred(:,year-min(year_range)+1) = gpinv(x0, a(1), a(2));

    t0 = rf_test.rf.(strcat('true_rf_', num2str(year)))(:);
    t1 = t0(t0>0);
    a = gpfit(t1);
    y_true(:,year-min(year_range)+1) = gpinv(x0, a(1), a(2));
    
    y_pred_stoc(:,year-min(year_range)+1) = rf_stoc.ypred_mean.(strcat('mean_', num2str(year)));
    
    p0 = rf_corr.rf_corrected.(strcat('pred_rf_corrected_', num2str(year)))(:);
    p1 = p0(p0>0);
    a = gpfit(p1);
    y_pred_corr(:,year-min(year_range)+1) = gpinv(x0, a(1), a(2));
end

figure()
plot(x,mean(y_true, 2),'g','LineWidth',3);hold on;
plot(x,mean(y_pred, 2),'r','LineWidth',3);hold on;
plot(x,mean(y_pred_stoc, 2),'b','LineWidth',3);hold on;
plot(x,mean(y_pred_corr, 2),'magenta','LineWidth',3);hold on;

L = mean(y_true, 2) - sqrt(mean((y_true - y_pred).^2,2));
H = mean(y_true, 2) + sqrt(mean((y_true - y_pred).^2,2));
plot(x,L,'-','LineWidth',1);
plot(x,H,'-','LineWidth',1);
fill([x(:); flipud(x(:))] ,[L(:);flipud(H(:))],'r','FaceAlpha',0.1);

L = mean(y_true, 2) - sqrt(mean((y_true - y_pred_stoc).^2,2));
H = mean(y_true, 2) + sqrt(mean((y_true - y_pred_stoc).^2,2));
plot(x,L,'-','LineWidth',1);
plot(x,H,'-','LineWidth',1);
fill([x(:); flipud(x(:))] ,[L(:);flipud(H(:))],'b','FaceAlpha',0.1);

L = mean(y_true, 2) - sqrt(mean((y_true - y_pred_corr).^2,2));
H = mean(y_true, 2) + sqrt(mean((y_true - y_pred_corr).^2,2));
plot(x,L,'-','LineWidth',1);
plot(x,H,'-','LineWidth',1);
fill([x(:); flipud(x(:))] ,[L(:);flipud(H(:))],'magenta','FaceAlpha',0.1);

% hold off;
% legend('Measured', 'Predicted (Deterministic)', 'Predicted (Stochastic)', 'Predicted (Bias Corrected)', 'Location', 'southeast');
legend('Observation (Daymet)', 'Prediction (Deterministic)', 'Prediction (Stochastic)', 'Prediction (Bias-Corrected)', 'Location', 'southeast');
set(gca,'XScale','log','FontSize',30);
xlabel('Return Period (years)', 'Fontsize', 30);
ylabel('Rainfall (mm/day)', 'Fontsize', 30);
% title('Predicted vs True Rainfall', 'Fontsize', 22);
grid on