%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%  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 '<name_d18O_file>.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 '<name_d18O_file>.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('<name_d18O_file>.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('<name_d18O_file>.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);