%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%  BIVALVE GROWTH MODEL %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% SCRIPT CREATED: January 2017 (E. Judd)
% LAST CHANGED: September 2017 (E. Judd)

% THIS IS THE SCRIPT TO RUN THE BIVALVE GROWTH MODEL PRESENTED IN:
    % Judd, E.J., Wilkinson, B.H., and Ivany, L.C. (2017). The life and
    % time of clams: Derivation of intra-annual growth rates from
    % high-resolution oxygen isotope profiles. Palaeogeography,
    % Palaeoclimatology, Palaeoecology.
    % https://doi.org/10.1016/j.palaeo.2017.09.034

% PRIOR TO RUNNING THIS SCRIPT, YOU WILL NEED TO DOWNLOAD THE SCEUA AND 
% CCEUA FUNCTIONS FROM THE MATLAB FILE EXCHANGE. THE FUNCTION CAN BE FOUND 
% AT:
    %https://www.mathworks.com/matlabcentral/fileexchange/7671-shuffled-complex-evolution--sce-ua--method
    
% YOU WILL NEED TO MODIFY THE SCEUA SCRIPT: 
    % In L72, replace 'functn' with 'GrowthModelFit'
% YOU WILL ALSO NEED TO MODIFY THE CCEUA SCRIPT:
    %In L43, L49, and L55, replace 'functn' with 'GrowthModelFit'
    
% This script is set up to run 'O1MA.csv', which consists
% of three years of synthetic data. It can be modified to run sequentially
% sampled isotope profiles from accretionary organisms.
% Regardless of the number of sampled years, input datasets should consist
% of 3 columns:
    %Column 1: Distance data (example data units are mm)
    %Column 2: Oxygen isotope values
    %Column 3: Year markers (1 = growth band, 0 = not a growth band)
% Refer to the 'O1MA.csv' for formatting of the file.

% This script utilizes the shuffled complex evolution iterative method
% (Duan et al., 1992), and calls the functions 'sceua.m' and 'cceua.m'
% developed by Q.Duan (9/2004) for more information on this method see:
    % Duan, Q., Sorooshian, S., Gupta, V., 1992. Effective and efficient
    % global optimization for conceptual rainfall-runoff models. Water
    % Resour. Res. 28, 1015?1031. doi:10.1029/91WR02985
    

% This script also calls the function 'GrowthModelFit.m' during the model
% run portion of the script and utilizes the function 'GrowthModel.m' to
% plot results.

% Please cite Judd et al. (2017) when using this code.
% Please address any questions to Emily Judd at ejjudd@syr.edu

%% STEP 1: READ IN DATA AND DEFINE VARIABLES
%Read in shell data
shelldata = csvread('O1MA.csv');
distdata = shelldata(:,1);
otopedata = shelldata(:,2);
yearmarker = find(shelldata(:,3)==1);
csvwrite('YEARMARKER.csv',yearmarker)
%Define variables
JD = linspace(1,365,365); per = 365;
%Set upper and lower bounds for Tsine and Gsine values
%Order of parameters is: [Tamp,Tpos,Gamp,Gpha,Gpos,Gskw,Tpha]
bl = [0,-4,0,0,-10,0,0]; %lower bound
bu = [5,2,10,365,10,100,365]; %upper bound
%Set initial Tsine and Gsine parameters - these values should be within
%the bounds set above
Tamp = 3; %Temp. (o-tope) Amplitude
Tpos = 0; %Temp (o-tope) Position
Tpha = 0;  %Temp (o-tope) Phase
Gamp = 5;  %Growth Amplitude
Gpha = 50; %Growth Phase
Gpos = 5;  %Growth Position
Gskw = 50; %Growth Skewness
x0 = [Tamp,Tpos,Gamp,Gpha,Gpos,Gskw,Tpha]; %Start parameters
%Set parameters for the shuffled complex evolution model
%for more details see the sceua function
iseed = 1; %the random seed number (for repetetive testing purpose)
iniflg = 1; %flag for initial parameter array (=1, included it in
%initial population; otherwise, not included)
ngs = 250; %number of complexes (sub-populations)
maxn = 1000*ngs*13; %maximum number of function evaluations allowed during optimization
kstop = 25; %maximum number of evolution loops before convergency
pcento = 0.001; %percentage change allowed in kstop loops before convergency
peps=0.001;


%% STEP 2: RUN MODEL

for k = 1:length(yearmarker)-1
    % extract a single year worth of distance measurements and isotope
    % values and write into a .csv file
    Dsamp = distdata(yearmarker(k):yearmarker(k+1)-1);
    Osamp = shelldata(yearmarker(k):yearmarker(k+1)-1,2);
    csvwrite('Osamp.csv',Osamp)
    csvwrite('Dsamp.csv',Dsamp)
    % define the function to iterate
    % function: GrowthModelFit
    % value to compare MSE: Osamp
    % start parameters: x0
    f0 = GrowthModelFit(Osamp,x0);
    [bestx,bestf] = sceua(x0,bl,bu,maxn,kstop,pcento,peps,ngs,iseed,iniflg);
    ExampleDataBestX(k,:) = bestx;
    ExampleDataBestF(k) = bestf;
end

% save solutions and MSE values for the model run
csvwrite('BESTX.csv',ExampleDataBestX);
csvwrite('BESTF.csv',ExampleDataBestF);

%% STEP 4: TEMPORALLY ALIGN FUNCTIONS WITH DATE OF MAX T AND PLOT RESULTS
shelldata = csvread('BESTX.csv');
yearmarker = csvread('YEARMARKER.csv');
shelldata = csvread('O1MA.csv');
distdata = shelldata(:,1);
otopedata = shelldata(:,2);
SW = -1; JD = linspace(1,365,365); per = 365;
bestx = csvread('BESTX.csv');
maxtemp = 182.5; %define the date of known or assumed max temp
for j = 1:length(yearmarker)-1
    Dsamp = distdata(yearmarker(j):yearmarker(j+1)-1);
    Osamp = otopedata(yearmarker(j):yearmarker(j+1)-1);
    csvwrite(strcat('Osamp_',num2str(j),'.csv'),Osamp)
    csvwrite(strcat('Dsamp_',num2str(j),'.csv'),Dsamp)
    Tamp = bestx(j,1);
    Tpha = bestx(j,7);
    Tpos = bestx(j,2);
    Gamp = bestx(j,3);
    Gpha = bestx(j,4);
    Gpos = bestx(j,5);
    Gskw = bestx(j,6);
    [Omod(:,j),Rmod(:,j),Dmod(:,j),Omodsamp(1:length(Dsamp),j),...
        Dmodsamp(1:length(Dsamp),j),JDmodsamp(1:length(Dsamp),j)] = ...
        GrowthModel(Tamp,Tpha,Tpos,Gamp,Gpha,Gpos,Gskw,Osamp,Dsamp);
    Tmod(:,j) = 18.03*1000./((Omod(:,j)-(0.97002*SW-29.98))+32.42)-273; %Calcite temp. equation of Kim and ONeill (1997)
    Tsamp = 18.03*1000./((Osamp-(0.97002*SW-29.98))+32.42)-273;
    csvwrite(strcat('Tsamp_',num2str(j),'.csv'),Tsamp)
    %Phase shift so that max temp occurs on the known or assumed date
    %of max temp
    JDshift(j) = maxtemp-mean(find(Tmod(:,j)==max(Tmod(:,j))));
    JDmodshift(1:length(JD),j) = JD(1:length(JD))+JDshift(j);
    JDmodsampshift(1:length(Dsamp),j) = JDmodsamp(1:length(Dsamp),j)+JDshift(j);
    %Adjust values so values fall within 1 and 365
    for i = 1:length(JD)
        if JDmodshift(i,j)<1
            JDmodshift(i,j) = JDmodshift(i,j)+365;
        elseif JDmodshift(i,j)>365
            JDmodshift(i,j) = JDmodshift(i,j)-365;
        end
    end
    for i = 1:length(Dsamp)
        if JDmodsampshift(i,j)<1
            JDmodsampshift(i,j) = JDmodsampshift(i,j)+365;
        elseif JDmodsampshift(i,j)>365
            JDmodsampshift(i,j) = JDmodsampshift(i,j)-365;
        end
    end
end
csvwrite('Omod.csv',Omod);
csvwrite('Rmod.csv',Rmod);
csvwrite('Dmod.csv',Dmod);
csvwrite('Omodsamp.csv',Omodsamp);
csvwrite('Dmodsamp.csv',Dmodsamp);
csvwrite('JDmodsamp.csv',JDmodsamp);
csvwrite('Tmod.csv',Tmod);
csvwrite('JDmodshift.csv',JDmodshift);
csvwrite('JDmodsampshift.csv',JDmodsampshift);

%% STEP 5: COLLATE OUTPUT INTO TABLES
% Create model output file
% Load Dmod data
Dmod = csvread('Dmod.csv');
Dmod2 = [Dmod;repmat(0,1,size(Dmod,2))];
Dmod3 = Dmod2(:);
% Load Omod data
Omod = csvread('Omod.csv');
Omod2 = [Omod;repmat(0,1,size(Omod,2))];
Omod3 = Omod2(:);
% Load Tmod data
Tmod = csvread('Tmod.csv');
Tmod2 = [Tmod;repmat(0,1,size(Tmod,2))];
Tmod3 = Tmod2(:);
% Load Rmod data
Rmod = csvread('Rmod.csv');
Rmod2 = [Rmod;repmat(0,1,size(Rmod,2))];
Rmod3 = Rmod2(:);
% Load JDmodshift data
JDmodshift = csvread('JDmodshift.csv');
JDmodshift2 = [JDmodshift;repmat(0,1,size(JDmodshift,2))];
JDmodshift3 = JDmodshift2(:);

% Calculate cumulative shifted model and modelled shell age
JDmodcum = repmat([1:365,0],1,size(JDmodshift,2)) + 365 * repelem(0:[size(JDmodshift,2)-1],366);
JDmodcum((1:size(JDmodshift,2))*366) = 0;
JDmodshiftcum = JDmodshift;
for i = 1:size(JDmodshiftcum,2)
    if i > 1
        %Add a year in case the next year needs to be shifted (growth stop)
        JDmodshiftcum(:,i) = JDmodshiftcum(:,i) - 365*floor((JDmodshiftcum(1,i)-JDmodshiftcum(end,(i-1)))/365);
    end

    N = find(JDmodshift([2:end 1],1)<JDmodshift(:,1))+1; %Find values going into new year
    JDmodshiftcum(N:end,i) = JDmodshiftcum(N:end,i) + 365; %Add 365 days to the values that are in the next year
end
JDmodshiftcum2 = [JDmodshiftcum;repmat(0,1,size(JDmodshift,2))];
JDmodshiftcum3 = JDmodshiftcum2(:);

% Combine into model output matrix
model_output = [Dmod3, Omod3, Tmod3, Rmod3, JDmodshift3, JDmodshiftcum3, JDmodcum(:)];

% Complete record by extrapolating the model in two directions
% Load original data
shelldata = csvread('O1MA.csv');

% Extrapolate model to start of record
firstyear = model_output(1:365,:); % Isolate the first year of the output matrix
if size(model_output,1) > 366
    firstyear(:,1) = firstyear(:,1) - (model_output(367,1) - model_output(1,1)); % Subtract the record length of one year
else
    firstyear(:,1) = firstyear(:,1) - (model_output(365,1) - model_output(1,1)); % Subtract the record length of one year, catch records with less than 2 years.
end
firstyear(:,6:7) = firstyear(:,6:7) - 365; % Subtract the number of days in one year
firstyear = [firstyear;repmat(0,1,size(firstyear,2))]; % Add zeroes to separate newly added lines
% Add extra data to start of record
if firstyear(1,1) < 0
    firstyear2 = firstyear(size(find(firstyear(:,1)<0),1):end,:);
    model_output2 = cat(1,firstyear2,model_output);
else
    model_output2 = cat(1,firstyear,model_output); % If adding a whole year does not reach the start of the record, max one more year is added
    firstyear2 = firstyear;
    firstyear2(:,1) = firstyear2(:,1) - (model_output(367,1) - model_output(1,1)); % Subtract the record length of one year
    firstyear2(:,6:7) = firstyear2(:,6:7) - 365; % Subtract the number of days in one year
    firstyear2(end,:) = repmat(0,1,size(firstyear2,2)); % Refresh zeroes at end of record
    if fistyear2(1,1) < 0
        firstyear3 = firstyear2(size(find(firstyear2(:,1)<0),1):end,:);
        model_output2 = cat(1,firstyear3,model_output2);
    else
        model_output2 = cat(1,firstyear3,model_output2); % Add max one second year to extend the record, dont extrapolate further
    end
end

% Extrapolate model to end of record
lastyear = model_output((end-365):end,:); % Isolate the last year of the output matrix
if size(model_output,1) > 366
    lastyear(:,1) = lastyear(:,1) + (model_output((end-1),1) - model_output((end-367),1)); % Add the record length of one year
else
    lastyear(:,1) = lastyear(:,1) + (model_output((end-1),1) - model_output((end-365),1)); % Add the record length of one year, catch records with less than 2 years.
end
lastyear(:,6:7) = lastyear(:,6:7) + 365; % Add the number of days in one year
lastyear(end,:) = repmat(0,1,size(lastyear,2)); % Refresh zeroes at end of record
if lastyear(end-1,1) > shelldata(end,1)
    lastyear2 = lastyear(1:(size(find(lastyear(:,1)<=shelldata(end,1)),1)+1),:); % Add lines up until one line after the last depth value in the original data
    model_output3 = cat(1,model_output2,lastyear2);
else
    model_output3 = cat(1,model_output2,lastyear); % If adding a whole year does not reach the end of the record, max one more year is added
    lastyear2 = lastyear;
    lastyear2(:,1) = lastyear2(:,1) + (model_output((end-1),1) - model_output((end-367),1)); % Add the record length of one year
    lastyear2(:,6:7) = lastyear2(:,6:7) + 365; % Add the number of days in one year
    lastyear2(end,:) = repmat(0,1,size(lastyear2,2)); % Refresh zeroes at end of record
    if lastyear2(end-1,1) > shelldata(end,1)
        lastyear3 = lastyear2(1:(size(find(lastyear2(:,1)<=shelldata(end,1)),1)+1),:); % Add lines up until one line after the last depth value in the original data
        model_output3 = cat(1,model_output3,lastyear3);
    else
        model_output3 = cat(1,model_output3,lastyear3); % Add max one second year to extend the record, dont extrapolate further
    end
end

model_export = array2table(model_output3,'VariableNames',{'Dmod', 'Omod', 'Tmod', 'Rmod', 'JD_mod_shift', 'JD_mod_shift_cum', 'JD_mod_cum'}); % Prepare model results for export
writetable(model_export,'model_output.csv'); % Export model results

% Create sample output file

SW = -1;
model_output4 = model_output3(any(model_output3,2),:); % Remove zero lines for linear interpolation
model_output4 = model_output4([diff(model_output4(:,1))~=0],:); % Remove duplicated x-values for linear interpolation
Dsamp = shelldata(:,1); % Copy Dsamp from original data
Osamp = shelldata(:,2); % Copy Osamp from original data
Tsamp = 18.03*1000./((Osamp-(0.97002*SW-29.98))+32.42)-273; % Calculate Tsamp from Osamp following Kim and ONeil, 1997 
Rsamp = interp1(model_output4(:,1),model_output4(:,4),Dsamp,'linear'); % find Rsamp by linearly interpolating Rmod
Dmodsamp = Dsamp;
JDmodsampshift = interp1(model_output4(:,1),model_output4(:,5),Dsamp,'linear'); % find JDmodsampshift by linearly interpolating JDmodshift
JDmodsampshiftcum = interp1(model_output4(:,1),model_output4(:,6),Dsamp,'linear'); % find JDmodsampshiftcum by linearly interpolating JDmodshiftcum
JDmodsampcum = interp1(model_output4(:,1),model_output4(:,7),Dsamp,'linear'); % find JDmodsampcum by linearly interpolating JDmodcum
JDmodsamp = mod(JDmodsampcum, 366); % find JDmodsamp by calculating the modulo of JDmodsampcum
Omodsamp = interp1(model_output4(:,1),model_output4(:,2),Dsamp,'linear'); % find Omodsamp by linearly interpolating Osamp
Tmodsamp = 18.03*1000./((Omodsamp-(0.97002*SW-29.98))+32.42)-273; % Calculate Tmodsamp from Omodsamp following Kim and ONeil, 1997
samp_output = [Dsamp, Osamp, Tsamp, Rsamp, Dmodsamp, JDmodsamp, JDmodsampshift, JDmodsampshiftcum, JDmodsampcum, Omodsamp, Tmodsamp]; % Combine sample interpolations into one matrix

% Insert zero rows between years
INS = find(diff(samp_output(:,6))<0);
samp_output2 = samp_output;
for i = 1:size(INS)
    samp_output2 = cat(1,samp_output2(1:(INS(i)+i-1),:),repmat(0,1,size(samp_output,2)),samp_output2((INS(i)+i):end,:));
end

sample_export = array2table(samp_output2,'VariableNames',{'Dsamp', 'Osamp', 'Tsamp', 'Rsamp', 'D_mod_samp' 'JD_mod_samp' 'JD_mod_samp_shift', 'JD_mod_samp_shift_cum', 'JD_mod_samp_cum' 'O_mod_samp' 'T_mod_samp'}); % Prepare sample results for export
writetable(sample_export,'sample_output.csv'); % Export sample results

% Create model parameters export
BESTF = csvread('BESTF.csv'); % Load BESTF data
BESTX = csvread('BESTX.csv'); % Load BESTX data
par_output = [transpose(1:size(BESTF,2)),transpose(BESTF),BESTX];
par_export = array2table(par_output,'VariableNames',{'YEAR' 'BESTF', 'T_amp', 'T_pos', 'T_phase', 'GR_amp' 'GR_pos' 'GR_phase', 'GR_skew'}); % Prepare parameter results for export
writetable(par_export,'parameter_output.csv'); % Export parameter results

% Combine into one Excel sheet
writetable(model_export, 'combined_results.xlsx', 'Sheet', 1);
writetable(sample_export, 'combined_results.xlsx', 'Sheet', 2);
writetable(par_export, 'combined_results.xlsx', 'Sheet', 3);