% Prepare the workspace
clear; close all; clc;

%Load data.
real = load('Validation.stats_real_proxy.txt');
pseudo = load('Validation.stats_pseudo_proxy.txt');
realind = load('Validation.stats_real_proxy_MEAN_REG.txt');
pseudoind = load('Validation.stats_pseudo_proxy_MEAN_REG.txt');

% Column structure in the CFR files is as follows.
% Col1: Calibration Grid RE
% Col2: Calibration Mean RE
% Col3: Calibration  r-squared
% Col4: Validation Mean Temperature
% Col5: Validation Mean Temperature (Fitted)
% Col6: Validation Grid RE
% Col7: Validation Mean RE
% Col8: Validation Mean CE
% Col9: Validation r-squared

% Column structure in the Index files is as follows.
% Col1: Calibration RE
% Col2: Calibration r-squared
% Col3: Validation Mean Temperature
% Col4: Validation Mean Temperature (Fitted)
% Col5: Validation RE
% Col6: Validation CE
% Col7: Validation r-squared

%Compute histograms for set number of bins
[nr,binr] = hist(real,450);
[np,binp] = hist(pseudo,2000);
[nrind,binrind] = hist(realind,450);
[npind,binpind] = hist(pseudoind,2000);

%Determine the maximum values for the normalized histograms
rmax = max(nr/length(real));
pmax = max(np/length(pseudo));
rindmax = max(nrind/length(realind));
pindmax = max(npind/length(pseudoind));

%Make figure 1
figure;
subplot(2,2,1)
%Below, as in all subsequent figures, the histograms for the red noise
%experiments is scaled to have the maximum value of the real proxy
%histograms.
bar(binp,(rmax(1)/pmax(1))*np(:,1)/length(pseudo),1,'facecolor','r','edgecolor','r')
hold on
bar(binr,nr(:,1)/length(real),1,'facecolor','k','edgecolor','k')
xlabel('RE')
ylabel('Normalized Frequency')
title('CFR Calibration Grid RE')
axis([-2 1 0 0.1])

subplot(2,2,3)
bar(binp,(rmax(3)/pmax(3))*np(:,3)/length(pseudo),1,'facecolor','r','edgecolor','r')
hold on
bar(binr,nr(:,3)/length(real),1,'facecolor','k','edgecolor','k')
xlabel('r^2')
ylabel('Normalized Frequency')
title('CFR Calibration Grid r^2')
axis([0 1 0 0.1])

subplot(2,2,2)
bar(binp,(rmax(6)/pmax(6))*np(:,6)/length(pseudo),1,'facecolor','r','edgecolor','r')
hold on
bar(binr,nr(:,6)/length(real),1,'facecolor','k','edgecolor','k')
xlabel('RE')
ylabel('Normalized Frequency')
title('CFR Validation Grid RE')
axis([-2 1 0 0.1])

subplot(2,2,4)
bar(binp,(rmax(9)/pmax(9)*np(:,9))/length(pseudo),1,'facecolor','r','edgecolor','r')
hold on
bar(binr,nr(:,9)/length(real),1,'facecolor','k','edgecolor','k')
xlabel('r^2')
ylabel('Normalized Frequency')
title('CFR Validation Grid r^2')
axis([0 1 0 0.1])

%Make figure 2
figure;
subplot(2,2,1)
bar(binp,(rmax(2)/pmax(2))*np(:,2)/length(pseudo),1,'facecolor','r','edgecolor','r')
hold on
bar(binr,nr(:,2)/length(real),1,'facecolor','k','edgecolor','k')
xlabel('RE')
ylabel('Normalized Frequency')
title('CFR Calibration Mean RE')
axis([-2 1 0 0.1])

subplot(2,2,3)
bar(binpind,(rindmax(1)/pindmax(1))*npind(:,1)/length(pseudoind),1,'facecolor','r','edgecolor','r')
hold on
bar(binrind,nrind(:,1)/length(realind),1,'facecolor','k','edgecolor','k')
xlabel('RE')
ylabel('Normalized Frequency')
title('Index Calibration RE')
axis([-2 1 0 0.1])

subplot(2,2,2)
bar(binp,(rmax(7)/pmax(7))*np(:,7)/length(pseudo),1,'facecolor','r','edgecolor','r')
hold on
bar(binr,nr(:,7)/length(real),1,'facecolor','k','edgecolor','k')
xlabel('RE')
ylabel('Normalized Frequency')
title('CFR Validation Mean RE')
axis([-2 1 0 0.1])

subplot(2,2,4)
bar(binpind,(rindmax(5)/pindmax(5)*npind(:,5))/length(pseudoind),1,'facecolor','r','edgecolor','r')
hold on
bar(binrind,nrind(:,5)/length(realind),1,'facecolor','k','edgecolor','k')
xlabel('RE')
ylabel('Normalized Frequency')
title('Index Validation RE')
axis([-2 1 0 0.1])

%Make figure 3
figure;
subplot(2,2,1)
bar(binp,(rmax(8)/pmax(8))*np(:,8)/length(pseudo),1,'facecolor','r','edgecolor','r')
hold on
bar(binr,nr(:,8)/length(real),1,'facecolor','k','edgecolor','k')
xlabel('CE')
ylabel('Normalized Frequency')
title('CFR Validation Mean CE')
axis([-2 1 0 0.1])

subplot(2,2,2)
bar(binp,(rmax(5)/pmax(5))*np(:,5)/length(pseudo),1,'facecolor','r','edgecolor','r')
hold on
bar(binr,nr(:,5)/length(real),1,'facecolor','k','edgecolor','k')
y = 0:0.01:0.1;
x = -0.35*ones(1,11);
plot(x,y,'k--')
xlabel('Temperature Anomaly')
ylabel('Normalized Frequency')
title('CFR Validation Mean Temperature')
axis([-2 1 0 0.1])

subplot(2,2,3)
bar(binpind,(rindmax(6)/pindmax(6))*npind(:,6)/length(pseudoind),1,'facecolor','r','edgecolor','r')
hold on
bar(binrind,nrind(:,6)/length(realind),1,'facecolor','k','edgecolor','k')
xlabel('CE')
ylabel('Normalized Frequency')
title('Index Validation Mean CE')
axis([-2 1 0 0.1])

subplot(2,2,4)
bar(binpind,(rindmax(4)/pindmax(4)*npind(:,4))/length(pseudoind),1,'facecolor','r','edgecolor','r')
hold on
bar(binrind,nrind(:,4)/length(realind),1,'facecolor','k','edgecolor','k')
plot(x,y,'k--')
xlabel('Temperature Anomaly')
ylabel('Normalized Frequency')
title('Index Validation Mean Temperature')
axis([-2 1 0 0.1])

