%% smoothing splines for 1610 ice core data
% Error on splines generated through bootstrap-based monte carlo simulations of the splines
% written by m. osman (mo549@cam.ac.uk) motivated by code/data from a. king (amyking@bas.ac.uk)
% october 2023

% LD = Law Dome, 13 yr resolution, 3e-3 at 50yr freq.
% WD = WAIS Divide, 5 yr resolution, 1.3e-3 at 50yr freq.
% ST = Skytrain, 8 yr resolution, 2e-3 at 50 yr freq.
% Comp2 = WAIS and Skytrain Data, 8e-5 at 100 yr freq.
% Comp3 = Law Dome, WAIS and Skytrain Data, 8e-5 at 100 yr freq.
% STCH4 = Skytrain CH4, 0.3 yr resolution, 8e-5 at 50 yr freq.

clear

%% Adjustable params

mc = 10000; % number of monte carlo realisations, suggested minimum == 10000 for publication, 1000 fine for sensitivity testing
yUnc = true; % option to include analytical uncertainty in the MC realisations, input true (more conservative errors) or false
rswr = true; % if true, does Random Sampling With Replacement, if false, then without replacement (this results in a fraction of measurements being removed)
splinePV = 3e-3; % changeable based on splinepv test, spline should be ~ half amplitude height (change based on the dataset you want to output and save the spline and errors for)

%test splinepv with generated function

t = 0:13:1000; %set to resolution of original CO2/CH4 data
a = 3;%arbritray
f = 0.02; %0.01 is 100 years, 0.02 is 50 years
y = a*cos(2*pi*f*t);
test = csaps(t,y,splinePV);
figure(1)
plot(t,y);
hold on
fnplt(test);
hold off
xlabel('Time')
ylabel('Amplitude')
title('Test Function')

%% load in data and set paramters for csaps function and monte carlo

% law dome
LD.dat = table2array(readtable('LD_data.csv'));
LD.int = [-1900:1:0]'; %-450:1:-275
LD.n = size(LD.dat,1); % for convienience later on

% WAIS Divide
WD.dat = table2array(readtable('WD_data.csv'));
WD.int = [-1200:1:-50]'; %-450:1:-275
WD.n = size(WD.dat,1); % for convienience later on

% Skytrain
ST.dat = table2array(readtable('ST_data_V2.csv'));
ST.int = [-495:1:-265]';
ST.n = size(ST.dat,1); % for convienience later on

% Skytrain CH4
STCH4.dat = table2array(readtable('ST_CH4_data.csv'));
STCH4.int = [-520:1:-335]';
STCH4.n = size(STCH4.dat,1); % for convienience later on

% WAIS and Skytrain composite
Comp2.dat = table2array(readtable('WD_ST_data_reset.csv'));
Comp2.int = [-1200:1:-50]';
Comp2.n = size(Comp2.dat,1); % for convienience later on

% Law Dome, WAIS and Skytrain composite
Comp3.dat = table2array(readtable('LD_WD_ST_data_reset.csv'));
Comp3.int = [-1900:1:0]';
Comp3.n = size(Comp3.dat,1); % for convienience later on

%% Run bootstrap
% random sampling with replacement

% preallocate the output
LD.spline = nan(length(LD.int),mc); 
WD.spline = nan(length(WD.int),mc); 
ST.spline = nan(length(ST.int),mc); 
STCH4.spline = nan(length(STCH4.int),mc); 
Comp2.spline = nan(length(Comp2.int),mc); 
Comp3.spline = nan(length(Comp3.int),mc); 

h = waitbar(0,'Conducting Monte Carlo simulations ... please wait!'); % include a waitbar to gauge level of impatience in real time :)
for i = 1:mc
    
    % grab current set of indices -- I are the bootstrap indices (i.e., randomly sampled  
    % data indices with replacement), meaning data indices can repeat, i'll
    % save the proportation of unique data points ("prUn") in each bootstap realisation for your curiosity 
    % LD
    LD.I = sort( randi(LD.n, [LD.n, 1]), 'ascend'); LD.Iu = unique(LD.I); 
    LD.prUn(i,1) = length(LD.Iu) / length(LD.I) * 100; 
    LD.dati = LD.dat(LD.I,:); 
    % WD
    WD.I = sort( randi(WD.n, [WD.n, 1]), 'ascend'); WD.Iu = unique(WD.I); 
    WD.prUn(i,1) = length(WD.Iu) / length(WD.I) * 100; 
    WD.dati = WD.dat(WD.I,:); 
    % ST
    ST.I = sort( randi(ST.n, [ST.n, 1]), 'ascend'); ST.Iu = unique(ST.I); 
    ST.prUn(i,1) = length(ST.Iu) / length(ST.I) * 100; 
    ST.dati = ST.dat(ST.I,:); 
    % STCH4
    STCH4.I = sort( randi(STCH4.n, [STCH4.n, 1]), 'ascend'); STCH4.Iu = unique(STCH4.I); 
    STCH4.prUn(i,1) = length(STCH4.Iu) / length(STCH4.I) * 100; 
    STCH4.dati = STCH4.dat(STCH4.I,:); 
    % Comp2
    Comp2.I = sort( randi(Comp2.n, [Comp2.n, 1]), 'ascend'); Comp2.Iu = unique(Comp2.I); 
    Comp2.prUn(i,1) = length(Comp2.Iu) / length(Comp2.I) * 100; 
    Comp2.dati = Comp2.dat(Comp2.I,:); 
    % Comp3
    Comp3.I = sort( randi(Comp3.n, [Comp3.n, 1]), 'ascend'); Comp3.Iu = unique(Comp3.I); 
    Comp3.prUn(i,1) = length(Comp3.Iu) / length(Comp3.I) * 100; 
    Comp3.dati = Comp3.dat(Comp3.I,:); 

    % option to do random sampling without replacement -- this removes all
    % non-unique data points, i.e., those whose indices have been listed
    % more than once.  This truncates the number of datapoints to <= n ...
    if ~rswr % if false...
        LD.dati = LD.dati(LD.Iu,:); 
        WD.dati = WD.dati(WD.Iu,:); 
        ST.dati = ST.dati(ST.Iu,:); 
        STCH4.dati = STCH4.dati(STCH4.Iu,:); 
        Comp2.dati = Comp2.dati(Comp2.Iu,:); 
        Comp3.dati = Comp3.dati(Comp3.Iu,:); 
    end

    % option to include analytical uncertainty -- this randomly samples the y-uncertainty
    % for each data point, assuming the 4th column in .dat represents N ~ (0,1sigma)
    if yUnc 
        LD.dati(:,3) = normrnd(LD.dati(:,3),LD.dati(:,4),[size(LD.dati,1),1]);
        WD.dati(:,3) = normrnd(WD.dati(:,3),WD.dati(:,4),[size(WD.dati,1),1]);
        ST.dati(:,3) = normrnd(ST.dati(:,3),ST.dati(:,4),[size(ST.dati,1),1]);
        STCH4.dati(:,3) = normrnd(STCH4.dati(:,3),STCH4.dati(:,4),[size(STCH4.dati,1),1]);
        Comp2.dati(:,3) = normrnd(Comp2.dati(:,3),Comp2.dati(:,4),[size(Comp2.dati,1),1]);
        Comp3.dati(:,3) = normrnd(Comp3.dati(:,3),Comp3.dati(:,4),[size(Comp3.dati,1),1]);
    end

    % finally, run the spline function
    LD.spline(:,i) = csaps(-LD.dati(:,2),LD.dati(:,3),splinePV, LD.int);
    WD.spline(:,i) = csaps(-WD.dati(:,2),WD.dati(:,3),splinePV, WD.int);
    ST.spline(:,i) = csaps(-ST.dati(:,2),ST.dati(:,3),splinePV, ST.int);
    STCH4.spline(:,i) = csaps(-STCH4.dati(:,2),STCH4.dati(:,3),splinePV, STCH4.int);
    Comp2.spline(:,i) = csaps(-Comp2.dati(:,2),Comp2.dati(:,3),splinePV, Comp2.int);
    Comp3.spline(:,i) = csaps(-Comp3.dati(:,2),Comp3.dati(:,3),splinePV, Comp3.int);
    
    % update the waitbar
    waitbar(i / mc)

end
close(h);

%% Some quick plots! 
% option to output spline and +/- 1 sigma confidence intervals for each
% dataset *remember to have set the appropriate splinePV for the dataset

xRange = [200 500]; yRange = [270 290]; % plotting ranges (except CH4)
col1 = [0.9961    0.8980    0.8510]; % light red
col2 = [0.9882    0.5725    0.4471]; % medium red
col3 = [0.8706    0.1765    0.1490]; % dark red

h = figure; hold on;
set(h,'units','centimeters','position',[1.8,0.90,14,28]);
    set(h,'PaperPositionMode','auto','PaperOrientation','landscape','Color','W');
    ax = gca; ax.Visible = 'off';

% LD
subplot(6,1,1); 
    set(gca,'Color','none','xdir','reverse'); hold on
    xlim(xRange); ylim(yRange); 
    xlabel('Year BP'); ylabel('Law Dome CO_{2} (ppm)')
    % 95% CI 
    time = -LD.int; 
    %low = prctile(LD.spline,2.5,2); upp = prctile(LD.spline,97.5,2); 
    %h1 = fill([time(1); time; flipud([time; time(end)])],[upp(1); low; flipud([upp; low(end)])], col1); 
    % +/- 1sigma
    low = nanmean(LD.spline,2) - nanstd(LD.spline,0,2); upp = nanmean(LD.spline,2) + nanstd(LD.spline,0,2);
    h2 = fill([time(1); time; flipud([time; time(end)])],[upp(1); low; flipud([upp; low(end)])], col2); 
    % mean
    plot(time,nanmean(LD.spline,2),'-','Linewidth',2,'Color',col3); 
    % data
    errorbar(LD.dat(:,2),LD.dat(:,3),2.*LD.dat(:,4), 'k.')
    scatter(LD.dat(:,2),LD.dat(:,3), 20, 'filled','MarkerFaceColor',0.1.*[1 1 1])

% LDspline_export = [time nanmean(LD.spline,2) low upp]
% filename = 'LDBootstrap.xlsx';
% writematrix(LDspline_export,filename);

% WD
subplot(6,1,2); 
    set(gca,'Color','none','xdir','reverse'); hold on
    xlim(xRange); ylim(yRange); 
    xlabel('Year BP'); ylabel('WAIS CO_{2} (ppm)')
    % 95% CI 
    time = -WD.int; 
    % low = prctile(WD.spline,2.5,2); upp = prctile(WD.spline,97.5,2); 
    % h1 = fill([time(1); time; flipud([time; time(end)])],[upp(1); low; flipud([upp; low(end)])], col1); 
    % +/- 1sigma
    low = nanmean(WD.spline,2) - nanstd(WD.spline,0,2); upp = nanmean(WD.spline,2) + nanstd(WD.spline,0,2);
    h2 = fill([time(1); time; flipud([time; time(end)])],[upp(1); low; flipud([upp; low(end)])], col2); 
    % mean
    plot(time,nanmean(WD.spline,2),'-','Linewidth',2,'Color',col3); 
    % data
    errorbar(WD.dat(:,2),WD.dat(:,3),2.*WD.dat(:,4), 'k.')
    scatter(WD.dat(:,2),WD.dat(:,3), 20, 'filled','MarkerFaceColor',0.1.*[1 1 1])

% WDspline_export = [time nanmean(WD.spline,2) low upp]
% filename = 'WDBootstrap.xlsx';
% writematrix(WDspline_export,filename);

% ST
subplot(6,1,3); 
    set(gca,'Color','none','xdir','reverse'); hold on
    xlim(xRange); ylim(yRange); 
    xlabel('Year BP'); ylabel('Skytrain CO_{2} (ppm)')
    % 95% CI 
    time = -ST.int; 
    % low = prctile(ST.spline,2.5,2); upp = prctile(ST.spline,97.5,2); 
    % h1 = fill([time(1); time; flipud([time; time(end)])],[upp(1); low; flipud([upp; low(end)])], col1); 
    % +/- 1sigma
    low = nanmean(ST.spline,2) - nanstd(ST.spline,0,2); upp = nanmean(ST.spline,2) + nanstd(ST.spline,0,2);
    h2 = fill([time(1); time; flipud([time; time(end)])],[upp(1); low; flipud([upp; low(end)])], col2); 
    % mean
    plot(time,nanmean(ST.spline,2),'-','Linewidth',2,'Color',col3); 
    % data
    errorbar(ST.dat(:,2),ST.dat(:,3),2.*ST.dat(:,4), 'k.')
    scatter(ST.dat(:,2),ST.dat(:,3), 20, 'filled','MarkerFaceColor',0.1.*[1 1 1])

% STspline_export = [time nanmean(ST.spline,2) low upp]
% filename = 'STBootstrap.xlsx';
% writematrix(STspline_export,filename);

% Comp2
subplot(6,1,4); 
    set(gca,'Color','none','xdir','reverse'); hold on
    xlim(xRange); ylim(yRange); 
    xlabel('Year BP'); ylabel('Comp2 CO_{2} (ppm)')
    % 95% CI 
    time = -Comp2.int; 
    % low = prctile(Comp2.spline,2.5,2); upp = prctile(Comp2.spline,97.5,2); 
    % h1 = fill([time(1); time; flipud([time; time(end)])],[upp(1); low; flipud([upp; low(end)])], col1); 
    % +/- 1sigma
    low = nanmean(Comp2.spline,2) - nanstd(Comp2.spline,0,2); upp = nanmean(Comp2.spline,2) + nanstd(Comp2.spline,0,2);
    h2 = fill([time(1); time; flipud([time; time(end)])],[upp(1); low; flipud([upp; low(end)])], col2); 
    % mean
    plot(time,nanmean(Comp2.spline,2),'-','Linewidth',2,'Color',col3); 
    % data
    errorbar(Comp2.dat(:,2),Comp2.dat(:,3),2.*Comp2.dat(:,4), 'k.')
    scatter(Comp2.dat(:,2),Comp2.dat(:,3), 20, 'filled','MarkerFaceColor',0.1.*[1 1 1])

% Comp2spline_export = [time nanmean(Comp2.spline,2) low upp]
% filename = 'Comp2Bootstrap.xlsx';
% writematrix(Comp2spline_export,filename);

% Comp3
subplot(6,1,5); 
    set(gca,'Color','none','xdir','reverse'); hold on
    xlim(xRange); ylim(yRange); 
    xlabel('Year BP'); ylabel('Comp3 CO_{2} (ppm)')
    % 95% CI 
    time = -Comp3.int; 
    % low = prctile(Comp3.spline,2.5,2); upp = prctile(Comp3.spline,97.5,2); 
    % h1 = fill([time(1); time; flipud([time; time(end)])],[upp(1); low; flipud([upp; low(end)])], col1); 
    % +/- 1sigma
    low = nanmean(Comp3.spline,2) - nanstd(Comp3.spline,0,2); upp = nanmean(Comp3.spline,2) + nanstd(Comp3.spline,0,2);
    h2 = fill([time(1); time; flipud([time; time(end)])],[upp(1); low; flipud([upp; low(end)])], col2); 
    % mean
    plot(time,nanmean(Comp3.spline,2),'-','Linewidth',2,'Color',col3); 
    % data
    errorbar(Comp3.dat(:,2),Comp3.dat(:,3),2.*Comp3.dat(:,4), 'k.')
    scatter(Comp3.dat(:,2),Comp3.dat(:,3), 20, 'filled','MarkerFaceColor',0.1.*[1 1 1])

% Comp3spline_export = [time nanmean(Comp3.spline,2) low upp]
% filename = 'Comp3Bootstrap.xlsx';
% writematrix(Comp3spline_export,filename);

% STCH4
subplot(6,1,6); 
    set(gca,'Color','none','xdir','reverse'); hold on
    xlim(xRange); ylim([665 735]); 
    xlabel('Year BP'); ylabel('Skytrain CH_{4} (ppb)')
    % 95% CI 
    time = -STCH4.int; 
    % low = prctile(STCH4.spline,2.5,2); upp = prctile(STCH4.spline,97.5,2); 
    % h1 = fill([time(1); time; flipud([time; time(end)])],[upp(1); low; flipud([upp; low(end)])], col1); 
    % +/- 1sigma
    low = nanmean(STCH4.spline,2) - nanstd(STCH4.spline,0,2); upp = nanmean(STCH4.spline,2) + nanstd(STCH4.spline,0,2);
    h2 = fill([time(1); time; flipud([time; time(end)])],[upp(1); low; flipud([upp; low(end)])], col2); 
    % mean
    plot(time,nanmean(STCH4.spline,2),'-','Linewidth',2,'Color',col3); 
    % data
    e = errorbar(STCH4.dat(:,2),STCH4.dat(:,3),2.*STCH4.dat(:,4), 'k.'); e.CapSize = 0;
    scatter(STCH4.dat(:,2),STCH4.dat(:,3), 20, 'filled','MarkerFaceColor',0.1.*[1 1 1])

% STCH4spline_export = [time nanmean(STCH4.spline,2) low upp]
% filename = 'STCH4Bootstrap.xlsx';
% writematrix(STCH4spline_export,filename);

