
%% This code inverts fossil SI values to estimate paleo-CO2 using functions
% derived perviously using FDA with training data sets
% 
clear all

%% 
% Load the data with contains the 'fine' output from the
% simple_mono_smoother code. CO2fine and SI_outfine

load simple_out_g_030409.mat

% Load the uncertain fossil obs to invert.
% It is expecting a 5 column matrix
% column (1) is the number of leaves etc. for that obs.
% column (2) is mean of the observation to be inverted (fossil leaf SI)
% column (3) is the standard deviation of this observation
% column (4) is the age of the obs
% column (5) is any previous estimate of the inversion (i.e. paleo-CO2)
% Any column length (within reason!) can be used and is determined by the
% code
% Only columns 2 and 3 are used here and is required, the other information
% is used in later plotting routines.

foss_obs=csvread('ginkgo_foss_data_22.csv');

[numobs nn] = size(foss_obs); % determine number of obs

mean_std(numobs,6)=zeros; % prepare array

%%
% PDFs are plotted of the observations inverted (i.e. fossil leaf SI) these
% indicate whether there has been sufficient sampling based on the mean and
% stdev and this information has been used correctly. 
% PDFS are also plotted of the inversions (i.e. paleo-CO2)

figure;
[aa]=gcf;
figure;
[bb]=gcf;

ylab='SI (%)';
xlab='CO_2 (ppmv)';

%% The main DO loop

for iobs = 1:numobs

%%
% Generate the estimates of obs for inversion.
% Input the number of estimates of obs for inversion, I found in the SI
% case that 500 produced repeatable results.
zzsize=500;
obs=foss_obs(iobs,2);
obsconf=foss_obs(iobs,3);  
zz = normrnd(0,1,zzsize,1);
% Add random gaussian error to each estimate
UnObs = obs+(obsconf.*zz); 

[f,ix]=ksdensity(UnObs); % Calcualte PDF of the obs estimates used in the inversion

figure(aa);
subplot(6,5,iobs)
plot(ix,f,'b');hold on

%% The inversion
% Find the CO2 values from the 2000 FDA curves that correspond to each SI
% esitmate rounded to 4 decimal places. You might want to tweak the number
% of decimal places to use.

for i = 1
    found=CO2fine(find((UnObs(i))-0.0001<SI_outfine & SI_outfine <(UnObs(i))+0.0001 ));
    [nfind nf]=size(found);
    nfind;
    foundc=found;
    [nfindc nf]=size(foundc);
    nfindc;
end  

for i = 2:zzsize
    found=CO2fine(find((UnObs(i))-0.0001<SI_outfine & SI_outfine <(UnObs(i))+0.0001 ));
    [nfind nf]=size(found);
    nfind;
    foundc=cat(1,found,foundc);
    [nfindc nf]=size(foundc);
    nfindc;
end    


found_out(iobs,1:nfindc)=foundc;
found_count(iobs)=nfindc;

[f,ix]=ksdensity(foundc); % Calcualte PDF of found CO2
figure(bb);
subplot(6,5,iobs)
plot(ix,f,'b');hold on

%%
% Determine the output statistics for the found CO2 values. Note here we
% determine the 95% confidence interval so the output columns are
% column (1) is the mean
% column (2) is the standard deviation
% column (3) is the 2.5 percentile
% column (4) is the median
% column (5) is the 97.5 percentile
% column (6) is n, the number of CO2 values found to match each uncertain
% fossil leaf SI
% We also save the information about the PDFs so these can be reproduced
% without having to recalculate them

mean_std(iobs,1)=mean(foundc);
mean_std(iobs,2)=std(foundc);
found_quart = quantile(foundc,[0.025,0.5,0.975]);
mean_std(iobs,3)=found_quart(1);
mean_std(iobs,4)=found_quart(2);
mean_std(iobs,5)=found_quart(3);
mean_std(iobs,6)=found_count(iobs);
kernel_1(iobs,:)=ix;
kernel_2(iobs,:)=f;
clear found
end
 
%%
% Save the statistics for further use, including by plotting routines

save('CO2_SI_statistics_060409','mean_std','found_out','found_count','kernel_1','kernel_2');