% This routine performs the wiggle-matching algorithm between pairs of Δ14C
% datasets and estimates a likely timescale offset and its uncertainty.
% Please cite the original creators of this algorithm [1], its first application to ice cores [2]
% and the paper accompanying this code [3]. 
% This code was written by Giulia Sinnl and released 10/05/2022. Please, contact
% Giulia Sinnl or Florian Adolphi for any more information.
% The Δ14C datasets used in this code are to be cited if used: NorthGRIP [3], WDC [3], 
% and Hulu-Cave speleothem data [4,5].
% All ages are in years b2k (before 2000 CE).
% References:
% [1] Bronk Ramsey, C. B., van der Plicht, J., & Weninger, B. ‘Wiggle matching’radiocarbon dates. Radiocarbon, 43(2A), find(h.BinEdges==1000)1-find(h.BinEdges==1000)9. doi: 10.1017/S003find(h.BinEdges==1000)222000find(h.BinEdges==1000)248 ,  2001
% [2] Adolphi, F., & Muscheler, R. Synchronizing the Greenland ice core and radiocarbon timescales over the Holocene–Bayesian wiggle-matching of cosmogenic radionuclide records. Climate of the Past, 12(1), 15-30. https://doi.org/10.5194/cp-12-15-2016 ,  2016
% [3] Sinnl G., Adolphi F., Christl M., Welten K.C.,  Woodruff T.,  Caffee M.,  Svensson A.,  Muscheler R.,   Rasmussen S.O. (2022) Synchronizing ice-core and U/Th time scales in the Last Glacial Maximum using Hulu Cave 14C and  new 10Be measurements from Greenland and Antarctica  Clim. Past. [submitted Aug 2022]
% [4] Southon, J., Noronha, A. L., Cheng, H., Edwards, R. L., & Wang, Y. A high-resolution record of atmospheric 14C based on Hulu Cave speleothem H82. Quaternary Science Reviews, 33, 32-41. https://doi.org/10.1016/j.quascirev.2011.11.022 ,  2012
% [5] Cheng, H., Edwards, R. L., Southon, J., Matsumoto, K., Feinberg, J. M., Sinha, A., ... & Ning, Y. Atmospheric 14C/12C changes during the last glacial period from Hulu Cave. science, 362(6420), 1293-1297. DOI: 10.1126/science.aau0747 ,  2018
%
% INPUT variables:

load('Delta14C_modelled_detrended'); % modelled and detrended D14C data from ice core 10Be (concentrations and fluxes)
load('Delta14C_Hulu_measured_detrended'); % detrended D14C data from Hulu Cave (ages converted yrs b2k)

T_min=2e4; T_max=2.5e4; % time boundaries [yrs b2k]
ts=[-1000:1:2000]; % t_offset vector to be scanned by the algorithm [yrs]
dW=50; % window overlap [yrs]
dT=1300; % window length [yrs]
W=[19300:dW:23800]; % time-window vector (will be limited by T_min and T_max later)
i_w=[1:length(W)-2]; % index of time-window vectors
Tw=W+dT/2; % mid-window time vector (for plotting)
t_grip=grip_age_shift(t_icecore); % GRIP ages according to the timescale shift proposed in [3]

% OUTPUT variables:
% P = P(W,t_offset) : 2 dim probability matrix of the wiggle-matching algorithm
% T_Mode{icecore} = t*_offset(W, icecore) : best time-dependent timescale offset for each ice core
% T_Mode_low68{icecore},T_Mode_high68{icecore} = lower and upper 1σ boundaries of T_Mode
% T_shift_gr = final average likely offset of Greenland datasets
% T_shift_ant = final average likely offset of Antarctica datasets
% err_shift_gr_low = Monte-Carlo based lower uncertainty bound of T_shift_gr
% err_shift_gr_high = Monte-Carlo based upper uncertainty bound of T_shift_gr
% err_shift_ant_low = Monte-Carlo based lower uncertainty bound of T_shift_ant
% err_shift_ant_high = Monte-Carlo based upper uncertainty bound of T_shift_ant

% Select either the MSD (i_hulu=1) or the H82 (i_hulu=2) dataset:
i_hulu=1;
if i_hulu==1
    t_hulu=t_hulu_cheng;
    y_hulu=y_hulu_cheng;
    e_hulu=e_hulu_cheng;
    et_hulu=et_hulu_cheng;
else
    t_hulu=t_hulu_southon;
    y_hulu=y_hulu_southon;
    e_hulu=e_hulu_southon;
    et_hulu=et_hulu_southon;
end
%% Wiggle matching
icecore_select=[1,2,3,4,8,11,12,13,14]; % corresponding dataset labels are stored in D14c_label
P={};
% wiggle matching routine (will take some time to perform)
P=wiggle_matching_routine(icecore_select,D14c_label,ts,W,dT,t_hulu,t_grip,t_icecore,D14c_modelled_detrended,D14c_error,y_hulu,e_hulu);

%% time dependent mode of P, with 68% confidence levels:
[T_Mode,T_Mode_low68,T_Mode_high68]=detect_mode_offset(P,D14c_label,icecore_select,i_w,t_grip,t_icecore,ts);

% Visually determined averaging boundaries for Greenland/Antarctica timescale offset:

T1_gr= 20900; T2_gr = 22000;
T1_ant= 20800; T2_ant = 21700;

%% Monte Carlo routine: Likely timescale offsets and errors

% Greenland
N_iter=10000;

mu_sample=cat(1,T_Mode(icecore_greenland)); % combine ice core T_Mode
sigma_low_sample=cat(1,T_Mode_low68(icecore_greenland)); % combine ice core T_Mode_low68
sigma_high_sample=cat(1,T_Mode_high68(icecore_greenland)); % combine ice core T_Mode_high68

sampling_idx=([1:length(mu_sample)]); % re-indexing ice cores
x= nan* zeros(length(Tw(i_w)),length(sampling_idx)*N_iter);

for i_Tw=find(Tw(i_w)>=T1_gr & Tw(i_w)<=T2_gr) % selecting wiggle-matching output within T1 and T2
    k=1;
    for ik=(sampling_idx)
        x_rand=[];
        j=1;
        % model P as a Gaussian with an average sigma:
        n_gauss=normpdf(ts,mu_sample{ik}(i_Tw),(sigma_low_sample{ik}(i_Tw)+sigma_high_sample{ik}(i_Tw))/2);
        n_gauss=n_gauss/max(n_gauss); % normalize between 0 and 1
        
        for iter=1:N_iter
            
            i_ts=randi(length(ts),1); % select a random offset index
            y_rand=rand(1,1); % select a random value between 0 and 1
            if y_rand<=n_gauss(i_ts) % if y<p(x)
                x_rand(j)=ts(i_ts); % create random variable, distributed according to n_gauss
                j=j+1;
            end
            
        end
        
        x(i_Tw, k+(0:length(x_rand)-1))=x_rand; % store them all in one big matrix
        k=k+1;
    end
end
hist_data_gr=x'; 

% Repeat for WDC

N_iter=10000;
mu_sample=cat(1,T_Mode(icecore_antarctica));
sigma_low_sample=cat(1,T_Mode_low68(icecore_antarctica));
sigma_high_sample=cat(1,T_Mode_high68(icecore_antarctica));

sampling_idx=([1:length(mu_sample)]);
x=nan*zeros(length(Tw(i_w)),length(sampling_idx)*N_iter);
for i_Tw=find(Tw(i_w)>=T1_gr & Tw(i_w)<=T2_gr)
    k=1;
    for ik=(sampling_idx)
        x_rand=[];j=1;
        n_gauss=normpdf(ts,mu_sample{ik}(i_Tw),(sigma_low_sample{ik}(i_Tw)+sigma_high_sample{ik}(i_Tw))/2);
        n_gauss=n_gauss/max(n_gauss);
        
        for iter=1:N_iter
            i_ts=randi(length(ts),1);
            y_rand=rand(1,1);
            if y_rand<=n_gauss(i_ts)
                x_rand(j)=ts(i_ts);j=j+1;
            end
        end
        x(i_Tw,k+(0:length(x_rand)-1))=x_rand;
        k=k+1;
    end
end
hist_data_ant=x';

%% Figures 

% colors for the ice cores:
color_y([1,5,6,7,8,9])={'k'};
color_y([2,10,11])={'r'};
color_y([3,12])={'b'};
color_y([4,13])={[0.016439210035888   0.591907149991785   0.568239306327787]};
color_y(14)={'m'};

%% Panel of T_Mode of Greenland ice cores

figure;
ax=subplot(2,2,1);hold on;
l=legend('orientation','horizontal','EdgeColor','w','NumColumns',3);
ax.InnerPosition(4)=0.7*ax.InnerPosition(4);
grid('on');

ax.Title.String='a) '; ax.Title.FontWeight='normal';
ax.XLabel.String='U/Th Age [ka b2k]';
ax.YLabel.String='t^*_{offset} [years]';

h=errorbar(21500,-800,[],[],-dT/2,+dT/2,'o','markersize',5,'color','k','markerfacecolor','b','parent',ax,'displayname','window width','linewidth',2);

xlim([20500 22500])
ylim([-1000,2000])

lw=1;
M_gr=[];
for i=[1,2,3,8,11,12,14]
    
    lw=1;
    if ismember(i,[8,11,12,14])
        lw=2;
    end
    s=plot(Tw(i_w),T_Mode{i},'-','color',color_y{i},'parent',ax,'displayname',label{i},'linewidth',lw);
    
end

% yellow patch
ylim1=get(gca,'YLim');
patch([T1_gr T2_gr T2_gr T1_gr],[ylim1(1) ylim1(1) ylim1(2) ylim1(2)],'y','handlevisibility','off','FaceAlpha',0.2,'EdgeColor','none')

ax.YAxis(1).TickValues=[-1000,0,125,400,675,1000,2000];
ax.XAxis.TickLabels=sprintfc('%.1f',(ax.XAxis.TickValues'/1000));

set(gcf, 'Units', 'centimeters','Position',[20 0 25 21]);
l.Position(2)=ax.Position(2)+ax.Position(4)+0.02;
l.Position(1)=ax.Position(1);
pos_ax_a=ax.Position;
%% Monte-Carlo histogram Greenland
ax=subplot(2,2,2);hold on;
ax.Position(4)=pos_ax_a(4);
ax.Title.String='a.1) ';ax.Title.FontWeight='normal';
ax.XLabel.String='t^*_{offset} [years]';
ax.YLabel.String='p(t*_{offset}) [arb. units]';

BW=50;
h=histogram(hist_data_gr(:,Tw(i_w)>=T1_gr & Tw(i_w)<=T2_gr),'BinWidth',BW, 'Normalization','probability','FaceColor','m','EdgeColor','k');

% the final mean is at the mode of the histogram:
[~,i_M]=max(h.Values);
T_shift_gr(i_hulu)=h.BinEdges(i_M)+BW/2;
xline(T_shift_gr(i_hulu),'k','linewidth',2);

Integral_right=h.Values(i_M)*BW/2+BW*sum(h.Values(i_M+1:end),'omitnan'); % integral on the right of the histogram
Integral_left=h.Values(i_M)*BW/2+BW*sum(h.Values(1:i_M-1),'omitnan'); % integral on the left of the histogram
Integral_tot=Integral_left+Integral_right;
if Integral_tot ~= sum(h.Values)*BW
    'something is wrong with the integral'
end

% calculate position of sigma+ and sigma- of the histogram:
index_sigma_left=0;
Integral_sigma_left=0;
while Integral_sigma_left<=(1-0.68)/2*Integral_tot % expand the left tail of histogram until it is 16%
    index_sigma_left=index_sigma_left+1;
    Integral_sigma_left=sum(BW*h.Values(1:index_sigma_left),'omitnan');
    
    if (index_sigma_left)==i_M
        break;
    end
end
if Integral_sigma_left>=(1-0.68)/2*Integral_tot
    'Integral_sigma_left is >= than 16%'
    index_sigma_left=index_sigma_left-1;
end
err_shift_gr_low(i_hulu)=h.BinEdges(index_sigma_left)+BW/2;
xline(err_shift_gr_low(i_hulu),'r','linewidth',1);

index_sigma_right=length(h.Values);
Integral_sigma_right=0;
while Integral_sigma_right<=(1-0.68)/2*Integral_tot % expand the right tail of histogram until it is 16%
    index_sigma_right=index_sigma_right-1;
    Integral_sigma_right=sum(BW*h.Values(index_sigma_right:end),'omitnan');
    
    if (index_sigma_right)==i_M
        break;
    end
end
if Integral_sigma_right>=(1-0.68)/2*Integral_tot
    'Integral_sigma_right is >= than 16%'
    index_sigma_right=index_sigma_right+1;
end
err_shift_gr_high(i_hulu)=h.BinEdges(index_sigma_right)-BW/2;
xline(err_shift_gr_high(i_hulu),'r','linewidth',1)

an=annotation('textbox','position',[0.1 0.9 0.1 0.1],'String',['Offset= ' num2str(T_shift_gr(i_hulu)) ' [' num2str(err_shift_gr_low(i_hulu)) '-' num2str(err_shift_gr_high(i_hulu)) '] years'])
pos_ax_a1=get(ax,'position');
an.Position=[pos_ax_a1(1) pos_ax_a1(2)+0.8*pos_ax_a1(4) pos_ax_a1(3) 0.2*pos_ax_a1(4)];
    
xlim([-1000,2000])
ylim([0,1.2*max(h.Values)])
%% Panel of T_Mode  Antarctica ice cores

ax2=subplot(2,2,3);hold on;
ax2.InnerPosition(4)=0.7*ax2.InnerPosition(4);

grid('on');
ylim([-1000 1000])
ax2.Title.String='b) ';ax2.Title.FontWeight='normal';
l=legend('orientation','horizontal','EdgeColor','w','NumColumns',3);
l.Position(1)=ax2.Position(1);
l.Position(2)=ax2.Position(2)+ax2.Position(4)+0.04;

ylim([-1000,1000]);

ax2.XLabel.String='U/Th Age [ka b2k]';
ax2.YLabel.String='t^*_{offset} [years]';

h=errorbar(21500,-800,[],[],-dT/2,+dT/2,'o','markersize',5,'color','k','markerfacecolor','b','parent',ax2,'displayname','window width','linewidth',2);

for i=[4,13]
    if i==4 | i==13
        lw=1;
        if i==4
            label{i}=' WDC conc.';lw=1;
        else
            label{i}='WDC flux';lw=2;
        end
s=plot(Tw(i_w),T_Mode{i},'-','color',color_y{i},'parent',ax2,'displayname',label{i},'linewidth',lw);
         
    end
end

ylim1=get(gca,'YLim');
patch([T1_ant T2_ant T2_ant T1_ant],[ylim1(1) ylim1(1) ylim1(2) ylim1(2)],'y','handlevisibility','off','FaceAlpha',0.2,'EdgeColor','none')

xlim([20500 22500])

ax2.YAxis(1).TickValues=[-1000,-125,0,225,575,1000];
ax2.XAxis.TickLabels=sprintfc('%.1f',(ax2.XAxis.TickValues'/1000));

pos_ax_b=ax2.Position;
%% Monte-Carlo histogram Antarctica
ax2=subplot(2,2,4);hold on;
ax2.Title.String='b.1) ';ax2.Title.FontWeight='normal';
ax2.XLabel.String='t^*_{offset} [years]';
ax2.YLabel.String='p(t*_{offset}) [arb. units]';
ax2.Position(4)=pos_ax_b(4);

BW=50;
h=histogram(hist_data_ant(:,Tw(i_w)>=T1_ant & Tw(i_w)<=T2_ant),'BinWidth',BW, 'Normalization','probability');
h.FaceColor=color_y{4};
h.EdgeColor='k';

[~,i_M]=max(h.Values);
T_shift_ant(i_hulu)=h.BinEdges(i_M)+BW/2;
xline(T_shift_ant(i_hulu),'k','linewidth',2);

Integral_right=h.Values(i_M)*BW/2+BW*sum(h.Values(i_M+1:end),'omitnan'); % integral on the right of the histogram
Integral_left=h.Values(i_M)*BW/2+BW*sum(h.Values(1:i_M-1),'omitnan'); % integral on the left of the histogram
Integral_tot=Integral_left+Integral_right;
if Integral_tot ~= sum(h.Values)*BW
    'something is wrong with the integral'
end

% calculate position of sigma+ and sigma-:
index_sigma_left=0;
Integral_sigma_left=0;
while Integral_sigma_left<=(1-0.68)/2*Integral_tot
    index_sigma_left=index_sigma_left+1;
    Integral_sigma_left=sum(BW*h.Values(1:index_sigma_left),'omitnan');
    
    if (index_sigma_left)==i_M
        break;
    end
end
if Integral_sigma_left>=(1-0.68)/2*Integral_tot
    'Integral_sigma_left is >= than 16%'
    index_sigma_left=index_sigma_left-1;
end
err_shift_ant_low(i_hulu)=h.BinEdges(index_sigma_left)+BW/2;
xline(err_shift_ant_low(i_hulu),'r','linewidth',1);

index_sigma_right=length(h.Values);
Integral_sigma_right=0;
while Integral_sigma_right<=(1-0.68)/2*Integral_tot
    index_sigma_right=index_sigma_right-1;
    Integral_sigma_right=sum(BW*h.Values(index_sigma_right:end),'omitnan');
    
    if (index_sigma_right)==i_M
        break;
    end
end
if Integral_sigma_right>=(1-0.68)/2*Integral_tot
    'Integral_sigma_right is >= than 16%'
    index_sigma_right=index_sigma_right+1;
end
err_shift_ant_high(i_hulu)=h.BinEdges(index_sigma_right)-BW/2;
xline(err_shift_ant_high(i_hulu),'r','linewidth',1);

an=annotation('textbox','position',[0.1 0.9 0.1 0.1],'String',['Offset= ' num2str(h.BinEdges(i_M)+BW/2) ' [' num2str(err_shift_ant_low(i_hulu)) '-' num2str(err_shift_ant_high(i_hulu)) '] years'])
pos_ax_b1=get(ax2,'position');
an.Position=[pos_ax_b1(1) pos_ax_b1(2)+0.8*pos_ax_b1(4) pos_ax_b1(3) 0.2*pos_ax_b1(4)];

xlim([-1000,2000]);
ylim([0 1.2*max(h.Values)]);

%% display results:
'Wiggle-matching algorithm terminated:'
disp(['Hulu-GICC05 offset = ' , num2str(T_shift_gr(i_hulu),3) ' years [', num2str([err_shift_gr_low,err_shift_gr_high]) '] (68% conf.)']);
disp(['Hulu-WD2014 offset = ' , num2str(T_shift_ant(i_hulu),3) ' years [', num2str([err_shift_ant_low,err_shift_ant_high]) '] (68% conf.)']);

%% figure of D14C data, shifted by the established offset

figure; hold on;
title('Figure 9a: Radionuclide data used for wiggle matching plotted on the synchronized time scales.')
l=legend();

h1=errorbar(t_hulu_cheng,y_hulu_cheng,e_hulu_cheng/2,'-o','markersize',2,'markerfacecolor',[0 10 49]/100,'color',[0 10 49]/100,'linewidth',1.5,'displayname','MSD (Hulu cave)');
h2=errorbar(t_hulu_southon,y_hulu_southon,e_hulu_southon/2,'-o','markersize',2,'markerfacecolor',[0 100 0]/100,'color',[0 100 0]/100,'linewidth',1.5,'displayname','H82 (Hulu cave)');
h1.CapSize=2; h2.CapSize=2;

plot(t_icecore+T_shift_gr(i_hulu),D14c_modelled_detrended{8},'DisplayName','NorthGRIP flux','color','k','linewidth',1);
plot(t_grip+T_shift_gr(i_hulu),D14c_modelled_detrended{11},'DisplayName','GRIP corrected flux','color','r','linewidth',1);
plot(t_icecore+T_shift_ant(i_hulu),D14c_modelled_detrended{13},'DisplayName','WDC flux','color','b','linewidth',1);
plot(t_icecore+T_shift_ant(i_hulu),D14c_modelled_detrended{14},'DisplayName','Flux stack','color','m','linewidth',1);
set(gcf, 'Units', 'centimeters','Position',[0 0 18 15]);

%% functions
function P=wiggle_matching_routine(i_icecore,label,ts,W,dT,t_hulu_res,t_grip,t_res,c14_d,errc14,y_hulu_res,e_hulu_res)
% iterate wiggle-matching for every selected ice-core dataset
for i=[i_icecore]
 f = waitbar(0,['Please be patient. Performing wiggle matching on ice core ' label{i} ' (' num2str(find(i_icecore==i)) '/' num2str(length(i_icecore)) ')']);
 
    if i==2 | i==10 | i==11
        t_icecore=t_grip;
    else
        t_icecore=t_res;
    end
    for iw=[1:length(W)]
        waitbar([iw/length(W)],f)
       
        T1=W(iw);
        T2=W(iw)+dT;        
        samplingrate= mean(diff(t_hulu_res));
% perform the actual wiggle matching on the selected data-window : 
            [P{iw,i},p1]=bayes_lag_prob_discont1(ts,T1,T2,t_icecore,c14_d{i}',...
                [t_icecore*0+errc14(i)],...% errc14 ‰ uncertainty for the ice core
                t_hulu_res,y_hulu_res,e_hulu_res,samplingrate,'plots',false);
    end
    delete(f)
end
end

function  [Mode,Mode_low68,Mode_high68]=detect_mode_offset(P,label,i_v,i_w,t_grip,t_res,ts)
for i=1:length(label)
    if ismember(i,[i_v])
        M_s{i}=[];k=1;
        if i==2 | i==10 | i==11
            t_icecore=t_grip;
        else
            t_icecore=t_res;
        end
        
        for iw=[i_w] % selected window indexes
            
            [~,i_M]=max(P{iw,i}); % detect the mode of P for every window W(iw)
            Mode{i}(k)=ts(i_M);
            
            % identify the main Probability lobe around the mode:
            [~,i_min_l]=min(P{iw,i}(1:i_M));
            [~,i_min_r]=min(P{iw,i}(i_M:end));
            i_min_r=i_min_r+i_M-1; %re-index to be the true index
            
            Mode_low68{i}(k)=0;
            Mode_high68{i}(k)=0;
            
            
            % detect the upper and lower FWHM boundaries of the probability
            i_left=i_M;
            while P{iw,i}(i_left)>=P{iw,i}(i_min_l)+0.5*(P{iw,i}(i_M)-P{iw,i}(i_min_l))
                i_left=i_left-1;
                if (i_left)==0
                    i_left=i_left+1;
                    break;
                end
                if i_left==i_M
                    break;
                end
            end
            i_left=i_left+1;
            
            i_right=i_M;
            while P{iw,i}(i_right)>=P{iw,i}(i_min_r)+0.5*(P{iw,i}(i_M)-P{iw,i}(i_min_r))
                i_right=i_right+1;
                if (i_right)==length(P{iw,i})+1
                    i_right=i_right-1;
                    break;
                end
                if (i_right)==i_min_r
                    break;
                end
            end
            i_right=i_right-1;

            Mode_low68{i}(k)=ts(i_M)-ts(i_left);
            Mode_high68{i}(k)=ts(i_right)-ts(i_M);
            
            if ts(i_M)==min(ts)
                Mode{i}(k)=nan;
            end
            k=k+1;
        end
    end
end
end
