%  ---------------------------------------------------------------------
%%       Program to fit curve between SI and CO2 concentrations
%       using monotone smoother of Ramsey (1998) taking into
%       account measurement errors in the training dataset.

%       This code relies completely upon that accompanying Ramsey 
%       and Silverman (2002) Applied Functional Data analysis: 
%       Methods and case studies, in particular Chapter 6.
%       Their MATLAB code is available from 
%       ftp://ego.psych.mcgill.ca/pub/ramsay/FDAfuns/Matlab
%       First download this and install using the installation
%       instructions there. 
%       You will have to edit the addpath below to reflect the location
%       you have installed this package.
%       
%       

%%  ---------------- input the data  ----------------------------------


% Add path to access FDA subroutines and functions. 

addpath ('../FDA') % Edit this to the location on your own system.

clear all;

% Load the data. The code is expecting 4 columns, 
% column (1) is mean of independent variable (CO2 conc.)
% column (2) is standard deviation of independent variable
% column (3) is mean of dependent variable (stomatal index)
% column (4) is standard deviation of dependent variable
% Any column length (within reason!) can be used and is determined by the code

a = csvread('ginkgo_data.csv');

CO2      = a(:,1);
CO2lo   = a(:,2);
SI      = a(:,3);
SIlo = a(:,4);
size_co2=size(CO2);
size_SI=size(SI);

%%
% Uncertainty is acoounted for by fitting curves to multiple pseudo-datasets
% generated using the mean and stdev of the actually observations.
% The number of psuedo-datasets to generate/curves to fit (nfit) needs to be 
% sufficant to generate repeatable curve statistics (i.e. mean and CI)
% For the SI example I found 2000 to be sufficant to do this

nfits=2000;

% Generate empty matricies to hold the pseudo-datasets

synobs_out=zeros(nfits,size_co2(1));
synobs_out2=zeros(nfits,size_co2(1));

%%
% The main DO loop - generate the pseudo-datasets and then fit curve

for i=1:nfits
    
% Generate new dataset by adding random guassian error taken from the
% uncertainty (standard deviation) in the observations. 

     z = normrnd(0,1,size_co2(1),1);
     newCO2 = CO2+(CO2lo.*z);
     z = normrnd(0,1,size_SI(1),1);
     newSI = SI+(SIlo.*z);
     synobs_out(i,:)=newSI;
     synobs_out2(i,:)=newCO2;

% Sort the observations in to accending order of independent variable
[n b]   = size(newCO2);
[nCO2,ind]  = sort(newCO2);
nSI  = newSI(ind);
rng  = [nCO2(1), nCO2(n)];

%%
% Set up the basis
% We use a B-spline of order 6 since we use a roughness penalty on the
% third derivative. 
% The number of knots to use will maybe require experimentation. Here I use
% a number equal to the number of observations. They are spaced at each data
% point location by passing the nCO2 values to create-bspline_basis.
% The number of basis functions this number of knots implies is n+6-2

nbasis   = n+6-2;
norder   = 6;
nCO2basis = create_bspline_basis(rng, nbasis, norder,nCO2);

% Wgt allows weighting of observations, in this case we have no premise for
% doing this so set to equal weighting. zmat is not used here, but could be
% used to allow variation in the coeffiecients over sampling points as a
% functino of some covariates

wgt  = ones(n,1);
zmat = wgt;

% Set up parameters for the monotone smooth
% lambda is the smoothing parameter determined through 12 fold cross
% validation

Lfd      = int2Lfd(2);
lambda   = 100;
nCO2fdPar = fdPar(nCO2basis, Lfd, lambda);

%  Carry out the monotone smooth

[Wfd, beta, Fstr, iternum, iterhist] = ...
    smooth_monotone(nCO2, nSI, nCO2fdPar);

%%
% When then evaluate the function for 5000 linearly spaced values of CO2
% between the min and max values - it is these 'curves' that we use in the
% following inversion step with the fossil observations.

CO2fine(i,:)  = linspace(nCO2(1),nCO2(n),5000)';
SI_out     = beta(1) + beta(2).*eval_mon(nCO2, Wfd);
SI_outfine(i,:) = beta(1) + beta(2).*eval_mon(CO2fine(i,:), Wfd);

end

%%
% The 'fine' curves are saved for use by the obs_inverter code
save('simple_out_g_030409','CO2fine','SI_outfine')
