%% Age-Depth Modelling in Complex Sedimentary Sequences
%
% Reference:
%
% Trauth, M.H., 2013, Age-Depth Modelling in Complex Sedimentary Sequences: The Problem of Extreme
% Variations of Sedimentation Rates and Hiatuses.
%
% Version 24 June 2013

%%
% Clear the workspace and load the data from the files. Add extra line at the end of |data| . 
clear, clc, close all
data = load('input_section.txt');
data(end+1,:) = [data(end,2) data(end,2) data(end,3)];
ages = load('input_ages.txt');

%%
% Define age |maxage| of the oldest slice of sediment in the section. This initial choice of the age
% can be modified iteratively during the experiment.
maxage = 401;

%% 
% The sediment column is devided into slices, each of them |slices| m thick. The algorithm aims at
% the determination of the most likely duration of the slices based on (1) the time difference
% between two age dates and (2) the unit deposition time of the sediment type.
slices = 0.1;

%%
% Use the |stairs| function to create a rectangular series of sediment types |stype| according to
% the base and top of the sediment layers in |data|. The section height or core depth are stored in
% the variable |depth|.
[depth,stype] = stairs(data(:,1),data(:,3),'o');

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(depth,stype,'r-',depth,stype,'bo')
xlabel('Height in Section [m]'), ylabel('Sediment Type 1, 2 or 3')
title('Sediment Types')

%%
% The interpolation of the sediment types |stype| on an evenly spaced grid |edepth| with the spacing
% |slices| requires the depth values in |depth| to be monotonically increasing. The |if| loop finds
% duplicate values in |depth| with different values in |stype|. The first |depth| value is reduced
% by 0.001 if |stype| is increasing. The first |depth| value is increased by 0.001 if |stype| is
% decreasing. Furthermore, subsequent |[depth,stype]| pairs are removed if subsequent |stype| data
% are identical since they are superfluous.
for i = 1 : length(depth)
    if i>=length(depth)-1, break, end
    if     stype(i) < stype(i+1)
           depth(i)   = depth(i)   -0.001;    
    elseif stype(i) > stype(i+1)
           depth(i+1) = depth(i+1) +0.001;             
    elseif stype(i) == stype(i+1) && ...
           stype(i+1) == stype(i+2)
           depth(i+1) = [];
           stype(i+1) = [];
    end
end

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(depth,stype,'r-',depth,stype,'bo')
xlabel('Height in Section [m]'), ylabel('Sediment Type 1, 2 or 3')
title('Sediment Types')

%%
% Interpolate |stype| on an evenly spaced grid |edepth|. The new depth vector |edepth| runs from the
% minimum to maximum value of the original depth vector |depth| with the spacing |slices|.
edepth = min(depth):slices:max(depth);
edepth = edepth';
estype = interp1(depth,stype,edepth,'linear');
edepth(1,:) = [];
estype(1,:) = [];

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(depth,stype,'o',edepth,estype,'r')
xlabel('Height in Section [m]'), ylabel('Sediment Type 1, 2 or 3')
title('Sediment Types')

%%
% Gamma model of the unit deposition time for each sediment type 1, 2 and 3, stored in |estype|. We
% use a Gamma model with a shape factor 2 and different scale factors 0.5, 0.8 and 2 for the three
% different sediment types. The scale factors can be iteratively adjusted in several runs of the 
% model. The modes |md| are calculated for the distributions, which can be easily replaced by any
% other measure of central tendency if required.
xa = -10 : 0.01 : 20; xa = xa';
xa = slices .* xa;
shape = [2 0.5; 2 0.8; 2 2];
shape(:,2) = slices .* shape(:,2);
for i = 1:length(shape)
    xt(:,i) = xa; yt(:,i) = ...
        gampdf(xt(:,i), shape(i,1), shape(i,2));
    yt(:,i) = yt(:,i) ./ sum(yt(:,i));
    md(i,1) = (shape(i,1)-1)*shape(i,2);
end

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(xt,yt)
xlabel('Time Accumulation Rate [kyr/m]'), ylabel('Frequency')
title('Typical Time Accumulation Rates')
legend(num2str(shape/slices)), legend('boxoff')

%%
% Create a vector |etimeaccsed| of the time accumulation per sediment slice with the thickness
% |slices| according to the sediment type in |estype|. |etimeaccsed| is based on the modes of the
% Gamma distributions of the sediment types in |estype|. The second array |etimedissed| is the set
% of Gamma distributions of each sediment slice according to the sediment type in |estype|. To
% reduce computation time, memory is first allocated for |etimeaccsed| and |etimedissed| before the
% |for| loop is started.
etimeaccsed = zeros(length(edepth),1);
etimedissed = zeros(length(edepth),length(xa));
for i = 1 : length(edepth)
    etimeaccsed(i,1) = slices .* md(uint8(estype(i)));
    etimedissed(i,:) = yt(:,uint8(estype(i)));
end

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(edepth,etimeaccsed,'r')
xlabel('Height in Section [m]'), ylabel('Accumulated Time per Slice [kyr]')
title('Time Accumulation Per Sediment Slice')

%%
% Calculate the mode, and the first, second (=median) and third quartile. The modes are stored in
% |etimedissedmodes|. The three quartiles are stored in |q1|, |q2| and |q3|.
for i = 1 : size(etimedissed,1)
    etimedissedmodes(i,1) = xa(etimedissed(i,:) == max(etimedissed(i,:)));
    xs = xa; xs = xs';    
    cs = cumsum(etimedissed(i,:)); cs = cs';
    xs(cs==0) = [];
    cs(cs==0) = [];
    xs(cs==max(cs)) = [];
    cs(cs==max(cs)) = [];  
    q1(i,:) = interp1(cs,xs,0.25);
    q2(i,:) = interp1(cs,xs,0.5);
    q3(i,:) = interp1(cs,xs,0.75);
end

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(edepth,etimedissedmodes,'b',edepth,q1,'r:',edepth,q2,'r',edepth,q3,'r:')
xlabel('Height in Section [m]'), ylabel('Accumulated Time per Slice [kyr]')
title('Time Accumulation Per Sediment Slice/Sedimentation Rates/Modes')
legend('Mode','1st Quartile','Median','3rd Quartile'), legend('boxoff')

%%
% Calculating the cumulative time accumulation |cummodes| of the section.
cummodes = cumsum(etimedissedmodes);

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(edepth,-(-maxage+cummodes),'b', ...
              edepth,-(-maxage+cummodes-q1/slices),'r-', ...
              edepth,-(-maxage+cummodes+q3/slices),'r-'); hold on
errorbar(ages(:,1),ages(:,2),ages(:,3),'+k');
xlabel('Height in Section [m]'), ylabel('Total Accumulated Time [kyr]')
title('Time Accumulation Per Sediment Slice/Sedimentation Rates/Modes')
legend('Mode','1st Quartile','3rd Quartile','Location','East')
legend('boxoff')

%%
% Calculate the difference |ageint| of subsequent age dates in |ages| using |diff|. Since the errors
% are Gaussian, we can calculate the errors of the age differences using the Gaussian law of error
% propagation.
clear ageint
ageint(:,1) = ages(:,1);
ageint(1:end-1,2) = abs(diff(ages(:,2)));
ageint(1:end-1,3) = sqrt(ages(1:end-1,3).^2 + circshift(ages(1:end-1,3),-1).^2);

%%
% Calculate the depth intervals |depthint| between the ages by substracting the midpoints of the
% dated layers from themselves, circular shifted by one value. Delete the last midpoint difference
% because this is the difference between the first and the last midpoint after circular shifting.
% Devide the intervals |depthint| by the thickness of the sediment slices |slices| to get the number
% of slices |depthslc| between the dated layers. Calculate the time of deposition and its error per
% slice by deviding the age interval |ageint| by the number of slices |depthslc|. Then, the time of
% deposition is stored in columns 2 and 3 for the array |ageint| where the initial age intervals and
% their errors were stored, i.e. they are overwritten by the new values.
depthint = abs(ages(:,1) - circshift(ages(:,1),-1));
depthint = depthint(1:end-1);
depthslc = depthint / slices;
ageint(1:end-1,2) = ageint(1:end-1,2)./depthslc;
ageint(1:end-1,3) = ageint(1:end-1,3)./(sqrt(2)*depthslc);

%%
% Use the |stairs| function to convert the age intervals and errors per slice into a rectangular
% representation. Herein, |adepth1| and |adepth2| are the depth vectors and |aageint| and |aageerr|
% are the age intervals and errors for each slice. The last value is deleted as it contains an
% |NaN| in both |aageint| and |aageerr|.
[adepth1,aageint] = stairs(ageint(:,1),ageint(:,2));
[adepth2,aageerr] = stairs(ageint(:,1),ageint(:,3));
adepth1(end) = [];
adepth2(end) = [];
aageint(end) = [];
aageerr(end) = [];

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(adepth1,aageint)
xlabel('Height in Section [m]'), ylabel('Accumulated Time per Slice [kyr]')
title('Time Accumulation Per Sediment Slice/Age Difference')

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(adepth2,aageerr)
xlabel('Height in Section [m]'), ylabel('Accumulated Time Error per Slice [kyr]')
title('Time Accumulation Error Per Sediment Slice/Age Difference')

%%
% The interpolation of the age differences and errors in |aageint| and |aageerr| on an evenly spaced
% grid |adepth1| and |adepth2| requires depth values in |adepth1| and |adepth2| to be monotonically
% increasing. The |if| loop finds duplicate values in |adepth1| and |adepth2| with different values
% in |aageint| and |aageerr|. The first |adepth1| and |adepth2| value is reduced by 0.001 if 
% |aageint| and |aageerr| is increasing. The first |adepth1| and |adepth2| value is increased by
% 0.001 if |aageint| and |aageerr| is decreasing. Furthermore, subsequent |[adepth1,aageint]| and
% pairs |[adepth2,aageerr]| are removed if subsequent |aageint| and |aageerr| data are identical
% since they are superfluous.
for i = 1 : length(adepth1)
    if i>=length(adepth1)-1, break, end
    if     aageint(i) < aageint(i+1)
           adepth1(i)   = adepth1(i)   -0.001;  
    elseif aageint(i) > aageint(i+1)
           adepth1(i+1) = adepth1(i+1) +0.001; 
    elseif aageint(i) == aageint(i+1) && aageint(i+1) == aageint(i+2)
           adepth1(i+1) = [];
           aageint(i+1) = [];
    end
end
for i = 1 : length(adepth2)
    if i>=length(adepth2)-1, break, end
    if     aageerr(i) < aageerr(i+1)
           adepth2(i)   = adepth2(i)   -0.001;  
    elseif aageerr(i) > aageerr(i+1)
           adepth2(i+1) = adepth2(i+1) +0.001; 
    elseif aageerr(i) == aageerr(i+1) && aageerr(i+1) == aageerr(i+2)
           adepth2(i+1) = [];
           aageerr(i+1) = [];
    end
end

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(adepth1,aageint,'o-')
xlabel('Height in Section [m]'), ylabel('Accumulated Time per Slice [kyr]')
title('Time Accumulation Per Sediment Slice/Age Difference')

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(adepth2,aageerr,'o-')
xlabel('Height in Section [m]'), ylabel('Accumulated Time Error per Slice [kyr]')
title('Time Accumulation Error Per Sediment Slice/Age Difference')

%% 
% Interpolate the |aageint| and |aageerr| on an evenly spaced grid |edepth|. We use the same evenly
% spaced depth vector |edepth| running from minimum to maximum value of the original depth vector
% |depth| with the spacing |slices|.
eageint = interp1(adepth1,aageint,edepth,'linear','extrap');
eageerr = interp1(adepth2,aageerr,edepth,'linear','extrap');

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(edepth,eageint,'r')
xlabel('Height in Section [m]'), ylabel('Accumulated Time per Slice [kyr]')
title('Time Accumulation Per Sediment Slice/Age Difference')

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(edepth,eageerr,'r')
xlabel('Height in Section [m]'), ylabel('Accumulated Time Error per Slice [kyr]')
title('Time Accumulation Error Per Sediment Slice/Age Difference')

%%
% Create a vector |etimeaccage| of the time accumulation per sediment slice with the thickness
% |slices| according to the age dates in |eageint|. |etimeaccage| is based on the modes of the age
% difference distribution. The second array |etimedisage| is the set of Gaussian distributions of
% each sediment slice according to the age dates in |eageint|. To reduce computation time, memory is
% first allocated for |etimeaccsed| and |etimedisage| before the |for| loop is started.
etimeaccage = zeros(length(edepth),1);
etimedisage = zeros(length(edepth),length(xa));
for i = 1 : length(edepth)
       etimeaccage(i,:) = eageint(i,:);
       etimedisage(i,:) = normpdf(xa,eageint(i,1),eageerr(i,1));
       etimedisage(i,:) = etimedisage(i,:) ./ nansum(etimedisage(i,:),2);
end

%%
% Calculate the mode, and the first, second (=median) and third quartile. The modes are stored in
% |etimedissedmodes|. The three quartiles are stored in |q1|, |q2| and |q3|.
for i = 1 : size(etimedisage,1)
    etimedissedmodes(i,1) = xa(find(etimedisage(i,:) == max(etimedisage(i,:))));
    %modes(i,1) = eageint(i,1);
    xs = xa; xs = xs';    
    cs = cumsum(etimedisage(i,:)); cs = cs';
    xs(cs==0) = [];
    cs(cs==0) = [];
    xs(cs==max(cs)) = [];
    cs(cs==max(cs)) = [];  
    q1(i,:) = interp1(cs,xs,0.25);
    q2(i,:) = interp1(cs,xs,0.5);
    q3(i,:) = interp1(cs,xs,0.75);
end

q1(q1<0) = 0;
q2(q2<0) = 0;
q3(q3<0) = 0;

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(edepth,etimedissedmodes,'b',edepth,q1,'r:',edepth,q2,'r',edepth,q3,'r:')
xlabel('Height in Section [m]'), ylabel('Accumulated Time per Slice [kyr]')
title('Time Accumulation Per Sediment Slice/Age Differences/Modes')
legend('Mode','1st Quartile','Median','3rd Quartile')
legend('boxoff')

%%
% Calculating the cumulative time accumulation |cummodes| of the section.
cummodes = cumsum(etimedissedmodes);

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(edepth,maxage-cummodes,'b', ...
              edepth,maxage-cummodes+q1/slices,'r-', ...
              edepth,maxage-cummodes-q3/slices,'r-'); hold on
errorbar(ages(:,1),ages(:,2),ages(:,3),'+k');
xlabel('Height in Section [m]'), ylabel('Total Accumulated Time [kyr]')
title('Time Accumulation Per Sediment Slice/Age Differences/Modes')
legend('Mode','1st Quartile','3rd Quartile','Location','East')
legend('boxoff')

%%
% Calculate the overlap |etimeoverlap| of the Gaussian distributions |etimedisage| of the time
% difference between two subsequent ages and the Gamma distributions |etimedissed| of the unit
% deposition time of the sediments by using the lower of the two curves for each |xa| value.
% Since both distributions are normalized to one, the overlap varies between 0 and 1, or the percent
% overlap between 0 and 100%. Segments of the section with low percent overlap between the two
% estimates of the deposition time are regarded as the most likely location of an hiatus.
xd = xa; xd = xd';
etimeoverlap = zeros(size(etimedissed));
for i = 1 : length(edepth)
    etimeoverlap(i,etimedisage(i,:)<etimedissed(i,:)) = ...
        etimedisage(i,etimedisage(i,:)<etimedissed(i,:));
    etimeoverlap(i,etimedissed(i,:)<=etimedisage(i,:)) = ...
        etimedissed(i,etimedissed(i,:)<=etimedisage(i,:));
    etimeoverlappercent(i) = max(cumsum(etimeoverlap(i,:)));
end

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(edepth,100*etimeoverlappercent,'r')
xlabel('Height in Section [m]'), ylabel('Percent Overlap [%]')
title('Percent Overlap between Age Estimates from Ages and Sediment Types')

%%
% Combining both estimates |etimedisage| and |etimedissed| of the deposition time by multiplying the
% distributions. The product of the two distributions makes the assumption that the distributions
% are independent and equally valid sources of information about the true deposition time. Merging
% the distributions to one yields a more reliable estimate of the deposition time. After merging the
% new distribution is again normalized to one.
etimedissall = etimedisage .* etimedissed;
for i = 1 : length(edepth)
    etimedissall(i,:) = etimedissall(i,:) ./ sum(etimedissall(i,:));
end

%%
% Calculate the mode, and the first, second (=median) and third quartile. The modes are stored in
% |etimedissallmodes|. The three quartiles are stored in |q1|, |q2| and |q3|.
for i = 1 : size(etimedissall,1)
    etimedissallmodes(i,1) = xa(etimedissall(i,:) == max(etimedissall(i,:)));
    xs = xa; xs = xs';    
    cs = cumsum(etimedissall(i,:)); cs = cs';
    xs(cs==0) = [];
    cs(cs==0) = [];
    xs(cs==max(cs)) = [];
    cs(cs==max(cs)) = [];  
    q1(i,:) = interp1(cs,xs,0.25);
    q2(i,:) = interp1(cs,xs,0.5);
    q3(i,:) = interp1(cs,xs,0.75);
end

%%
% Calculating the cumulative time accumulation |cummodes| of the section.
cummodes = cumsum(etimedissallmodes);

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(edepth,maxage-cummodes,'b', ...
              edepth,maxage-cummodes+q1/slices,'r-', ...
              edepth,maxage-cummodes-q3/slices,'r-'); hold on
errorbar(ages(:,1),ages(:,2),ages(:,3),'+k');
xlabel('Height in Section [m]'), ylabel('Total Accumulated Time [kyr]')
title('Time Accumulation Per Sediment Slice/Combined Sources of Evidence/Modes')
legend('Mode','1st Quartile','3rd Quartile','Location','East')
legend('boxoff')

%%
% Display age model together with percent overlap. Obviously the percent overlap is very low between
% 33.9 and 36 m height, which is the most likely location of an hiatus.
figure('Color',[1 1 1],'Position',[50 1000 600 300]);
[axes1,line1,line2] = plotyy(edepth,maxage-cummodes, ...
                        edepth,100*etimeoverlappercent); hold on
line3        = plot(edepth,maxage-cummodes,'b', ...
                        edepth,maxage-cummodes+q1/slices,'r-', ...
                        edepth,maxage-cummodes-q3/slices,'r-');
line4        = errorbar(ages(:,1),ages(:,2),ages(:,3),'+k');
set(get(axes1(1),'Title'),...
  'String','Time Accumulation Per Sediment Slice/Combined Sources of Evidence/Modes/% Overlap')
set(get(axes1(1),'XLabel'),...
  'String','Height in Section [m]')
set(get(axes1(1),'YLabel'),...
  'String','Total Accumulated Time [kyr]')
set(get(axes1(2),'YLabel'),...
  'String','Percent Overlap [%]')
set(axes1(1),'YColor',[1 0 0])
set(axes1(2),'YColor',[0 0.5 0])
set(line1,'Color',[0 0 0])
set(line2,'Color',[0 0.5 0])
set(line3(1),'Color',[0 0 1])
set(line3(2:3),'Color',[1 0 0])
set(line4,'Color',[0 0 0],'LineStyle',':')
grid

%%
% Having determined the most likely location of a hiatus between 33.9 and 36 m height, introduce a
% hiatus of 2 kyr length at 35 m.
cummodes_h = cummodes;
cummodes_h(edepth>35)  = cummodes_h(edepth>35) + 2;

figure('Color',[1 1 1],'Position',[50 1000 600 300]);
plot(edepth,maxage-cummodes_h,'b', ...
              edepth,maxage-cummodes_h+q1/slices,'r-', ...
              edepth,maxage-cummodes_h-q3/slices,'r-'); hold on
errorbar(ages(:,1),ages(:,2),ages(:,3),'+k');
xlabel('Height in Section [m]'), ylabel('Total Accumulated Time [kyr]')
title('Time Accumulation Per Sediment Slice/Combined Sources of Evidence and Hiatuses/Modes')
legend('Mode','1st Quartile','3rd Quartile','Location','East')
legend('boxoff')

%%
% Merge age depth data into one a single variable |agemodel|.
agemodel(:,1) = edepth;
agemodel(:,2) = -maxage+cummodes_h;