close all; clear all; clc;

load Sula_and_NORDOYAN_last;


% figure
% clf
% plot(data((data(:,7)>=datenum('2000','yyyy') & data(:,7)<=datenum('2008','yyyy')),7),...
%     data((data(:,7)>=datenum('2000','yyyy') & data(:,7)<=datenum('2008','yyyy')),8),'.-k')
% set(gca,'XTick',xpts,...
%     'XTickLabel',datestr(xpts,'yyyy'))
% xlabel('time')
% ylabel('Wind speed, m/s')
% % legend('observations')
% editplot2
% print(gcf, '-depsc','Wind_TS.eps');
% close
% 
% figure
% clf
% plot(data((data(:,7)>=datenum('2000','yyyy') & data(:,7)<=datenum('2008','yyyy')),7),...
%     data((data(:,7)>=datenum('2000','yyyy') & data(:,7)<=datenum('2008','yyyy')),10),'.-k')
% set(gca,'XTick',xpts,...
%     'XTickLabel',datestr(xpts,'yyyy'))
% xlabel('time')
% ylabel(['Hs (total sea), m'])
% % legend('observations')
% editplot2
% print(gcf, '-depsc','Hs_swell_TS.eps');
% close

% figure;
% clf;
% plot(X(:,9),Y(:,9),'.k')
% % legend('observations')
% xlabel('\xi')
% ylabel('\eta','Rotation',0)%'Position',[-0.095 0.5])
% editplot_paper

epsilon = @(eta,x) x(1)*exp(-x(3)*(eta-x(2)).^x(4));
star = @(level,x) x(2)+(1/x(3).*(log(x(1))-log(level))).^(1/x(4));


T = [2 5 10 20 50 100];
pr = 1 - 1./T;

k_memory = [1 2 4 24 48 72 96];

[maxX, yrs_X, Xmom(1), Xmom(2), design_mom_X] = MOM_Gumb(X, pr);
[maxY, yrs_Y, Ymom(1), Ymom(2), design_mom_Y] = MOM_Gumb(Y, pr);

r2_annual = corrcoef(maxX,maxY);
r2_annual = r2_annual(2,1);

kendall_annual = corr(maxX,maxY,'type','Kendall');
sherman_annual = corr(maxX,maxY,'type','Spearman');

% if (yrs_X == yrs_Y) && (yrs_X > 30)
%     [a_mle_X, b_mle_X, design_mle_X] = MLE_Gumb(X, pr);
%     [a_mle_Y, b_mle_Y, design_mle_Y] = MLE_Gumb(Y, pr);
% end

x_tmp = X(:)';
x_tmp(isnan(x_tmp)) = [];
Npoints = numel(x_tmp);
sigma_1 = std(x_tmp);

y_tmp = Y(:)';
y_tmp(isnan(y_tmp)) = [];
sigma_2 = std(y_tmp);

r2_all = corrcoef(x_tmp,y_tmp);
r2_all = r2_all(2,1);

% kendall_all = corr(x_tmp',y_tmp','type','Kendall');
% sherman_all = corr(x_tmp',y_tmp','type','Spearman');

ro = [r2_annual, r2_all];

clear x_tmp y_tmp r2_annual r2 r2_eta kendall_annual sherman_annual ...
    kendall_eta sherman_eta kendall_all sherman_all


% k_memory = [8 16];
% k_memory = [2 3 8 16 24 32];
% N_X = 500;
% N_Y = 300;
% N_X = numel(unique(X))+100;
% N_Y = numel(unique(Y))+100;
N_X = 500;
N_Y = 500;

%% ACER_1 on X
flagACER = 0;

[ bl_X, acer_hat_X,  CI_X] = uniACER( X, k_memory, N_X, flagACER);
close
%% ACER_1 on Y

[ bl_Y, acer_hat_Y,  CI_Y ] = uniACER( Y, k_memory, N_Y, flagACER);
close
%% Bivariate ACER
[~,ind] = max(log(acer_hat_X{end}));
eta_0_X = bl_X{end}(ind);
clear ind

[~,ind] = max(log(acer_hat_Y{end}));
eta_0_Y = bl_Y{end}(ind);
clear ind

h = waitbar(0,'Calculating ACER functions, please wait');
for kk = 1:numel(k_memory)
    
    %     [bar_lev{kk}, acer_hat{kk}, CI{kk}] = ...
    %         ACER(x, length(x)*ones(size(x, 1),1), k_memoryory(kk), N_barrier, conf_level, flag);
    [ACER_hat_mean{kk}, CI_BI{kk}, bar_lev_X{kk}, bar_lev_Y{kk}] ...
        = BiACER_uncutted( X, Y, eta_0_X, eta_0_Y, N_X, N_Y, k_memory(kk), flagACER);
    waitbar(kk/length(k_memory));
end
close(h);
clear kk h 

delta = 0.9;
eta_end_X = 38.32;
eta_end_Y = 43.32;
for kk = 1:numel(k_memory)

[ACER_hat_mean{kk}, CI_BI{kk}, bar_lev_X{kk}, bar_lev_Y{kk} ] = ...
     cutter_plot( ACER_hat_mean{kk}, CI_BI{kk}, bar_lev_X{kk}, bar_lev_Y{kk}, eta_0_X, eta_0_Y, eta_end_X, eta_end_Y );
end
clear kk eta_0_X eta_0_Y;

%% surfaces and contours

% red orange yellow green blue gray
color = {[1 0 0], [0.2 1 0.2], [1,0.65,0], [1 0.9 0],[0 0 1],[0.45 0.45 0.45],[0.63 0.13 0.94],[0.4 0.58 0.93],[1 0 0],[1 0 0]};
line = 1.5*ones(1,20);
%     v = -4 : 0.25 : -1.5;

% for kk = 1:numel(k_memory)
%     minn(kk) = min(min(log10(ACER_hat_mean{kk})));
%     maxx(kk) = max(max(log10(ACER_hat_mean{kk})));
% end
% % a = roundn(max(minn),-2);
% % b = roundn(min(maxx),-2)-0.5;
% a = max(minn);
% b = min(maxx);
% lev = b:-0.15:a;
% clear minn maxx a b
% for ii = 1:numel(lev)
%     
%     v = [lev(ii) lev(ii)];
%     
%     figure;
%     clf
%     for kk = 1:numel(k_memory)
%         [X_mesh,Y_mesh]=meshgrid(bar_lev_X{kk},bar_lev_Y{kk});
%         
%         [ACERContours{kk}, cont_handle{kk}] = contour(Y_mesh,X_mesh,log10(ACER_hat_mean{kk}),v,...
%             'Color',color{kk},'LineWidth',line(kk));
%         set(cont_handle{kk},'TextList', roundn(get(cont_handle{kk},'LevelList'),-2),...
%             'LevelList', roundn(get(cont_handle{kk},'LevelList'),-2))
%         hold on
%         leg{kk} = ['k=' num2str(k_memory(kk))];
%     end
%     txt = clabel(ACERContours{1},cont_handle{1},'LabelSpacing',1000);
%     set(txt,'BackgroundColor',[1 1 .6],'Rotation',0,...
%         'Edgecolor',[.7 .7 .7],'FontName', 'Times New Roman', 'FontSize', 18);
%     txt = clabel(ACERContours{numel(k_memory)},cont_handle{numel(k_memory)},'LabelSpacing',1000);
%     set(txt,'BackgroundColor',[1 1 .6],'Rotation',0,...
%         'Edgecolor',[.7 .7 .7],'FontName', 'Times New Roman', 'FontSize', 18);
%     legend(leg)
%     xlabel('\xi')
%     ylabel('\eta','Rotation',0)
%     text('Units','pixels','Position',[721,250],'String',['L = ' num2str(lev(ii),'%6.2f')],'BackgroundColor',[1 1 .6],...
%         'Edgecolor',[.7 .7 .7],'FontName', 'Times New Roman', 'FontSize', 18)
%     editplot_contour
% %     axis([1 11 14 27])
%     
%     saveas(gcf,['contour_' num2str(ii) '.fig']);
%     %     set(gcf,'PaperPositionMode','auto');
%     %     set(gcf,'InvertHardcopy','off')
%     %     print(gcf, '-depsc', [num2str(abs(100*lev)) '_contour.eps']);
%     %     clear h txt X Y
%     %
%     set(gcf,'PaperPositionMode','auto');
%     set(gcf,'PaperUnits','centimeters','PaperSize',[21 15]);
%     set(gcf,'InvertHardcopy','off')
%     print(gcf,'-dpdf', ['contour_' num2str(ii) '.pdf']);
%     close
% end

% saveas(gcf,[num2str(abs(100*lev)) '_contour.fig']);
% set(gcf,'PaperPositionMode','auto');
% set(gcf,'InvertHardcopy','off')
% print(gcf, '-depsc', [num2str(abs(100*lev)) '_contour.eps']);
clear h txt kk leg v lev ii X_mesh Y_mesh

figure;
clf
for kk = 1:numel(k_memory)
    [X_mesh,Y_mesh] = meshgrid(bar_lev_X{kk},bar_lev_Y{kk});
    
    mesh(Y_mesh,X_mesh,log10(ACER_hat_mean{kk}),'EdgeColor',color{kk},'Marker','.','MarkerSize',0.5);
    leg{kk} = ['k=' num2str(k_memory(kk))];
    hold on
end
view([110,7])
legend(leg)

xlabel('\xi')
ylabel('\eta')
h = zlabel('$\hat{\mathcal{E}}_{k}$','Rotation',0);
set(h,'Interpreter','LaTex')
set(gca,'ZTick',[-4 -3 -2 -1],'ZTickLabel',[-4 -3 -2 -1])
editplot_surface
% set(gca,'ZTick',[-4 -3 -2 -1],'ZTickLabel',[-4 -3 -2 -1])
% editplot_surface
% saveas(gcf,'surfaces.fig');
% set(gcf,'PaperPositionMode','auto');
% set(gcf,'InvertHardcopy','off')
% print(gcf, '-depsc', 'surfaces.eps');
clear X_mesh Y_mesh kk leg
%% Choosing the one
k_choise = 1; % 96 1
level = ( -log(pr)*yrs_X ) / (Npoints-k_choise+1);

ind = find(k_memory == k_choise);

% load Sula_and_NORDOYAN/allACERs;

bl_X = bl_X{ind};
acer_hat_X = acer_hat_X{ind};
CI_X = CI_X{ind};

bl_Y = bl_Y{ind};
acer_hat_Y = acer_hat_Y{ind};
CI_Y = CI_Y{ind};

bar_lev_X = bar_lev_X{ind};
bar_lev_Y = bar_lev_Y{ind};
ACER_hat_mean = ACER_hat_mean{ind};
CI_BI = CI_BI{ind};

clear ind
%% optimization!!!!
eta_1_X = 15; % 24 15
eta_1_Y = 15; % 27 15
delta = 0.9;
% optimization on X
[fin_sol_X, ~, sol_X, acerXtmp, blXtmp] = Opt(X, bl_X, acer_hat_X, CI_X, eta_1_X, delta);
X_ACER_star = star(level,fin_sol_X);

acer_fit_Xtmp = epsilon(blXtmp, fin_sol_X);
figure
clf
semilogy(blXtmp, acerXtmp, '*k','MarkerSize',2)
hold on
semilogy(blXtmp, acer_fit_Xtmp, 'k')
xlabel('\xi')
ylabel('ACER^x_2(\xi)')
legend('\epsilon_2(\xi)','\epsilon_2^{fit}(\xi)','CI^{+}','CI^{-}')
editplot2;
% print(gcf, '-depsc','acer_X_fit.eps');

% optimization on Y
[fin_sol_Y, ~, sol_Y, acerYtmp, blYtmp] = Opt(Y, bl_Y, acer_hat_Y, CI_Y, eta_1_Y, delta);
Y_ACER_star = star(level,fin_sol_Y);

acer_fit_Ytmp = epsilon(blYtmp, fin_sol_Y);
figure
clf
semilogy(blYtmp, acerYtmp, '*k','MarkerSize',2)
hold on
semilogy(blYtmp, acer_fit_Ytmp, 'k')
xlabel('\eta')
ylabel('ACER^y_2(\eta)')
legend('\epsilon_2(\eta)','\epsilon_2^{fit}(\eta)','CI^{+}','CI^{-}')
editplot2;
% print(gcf, '-depsc','acer_Y_fit.eps');

% optimization on Bivariate data
[ACER_hat_mean, CI_BI, bar_lev_X, bar_lev_Y ] = ...
    cutter( ACER_hat_mean, CI_BI, bar_lev_X, bar_lev_Y, eta_1_X, eta_1_Y, delta );

condX = ismember(bl_X,bar_lev_X);
acer_hat_X = acer_hat_X(condX);
bl_X = bl_X(condX);
clear condX

condY = ismember(bl_Y,bar_lev_Y);
acer_hat_Y = acer_hat_Y(condY);
bl_Y = bl_Y(condY);
clear condY

acer_fit_X = epsilon(bl_X, fin_sol_X);
acer_fit_Y = epsilon(bl_Y, fin_sol_Y);

clear bl_X bl_Y

%% Gumbel assymetric logistic model

[fin_sol, fmin, pos, sol, Weights] =...
    Optim_Bivar(ACER_hat_mean, CI_BI, acer_fit_X, acer_fit_Y, ro(2));

% Gumb_assym_log_ACER_hat = Gumb_assym_log(acer_hat_X,acer_hat_Y,fin_sol(1),fin_sol(2),fin_sol(3));
Gumb_assym_log_ACER_fit = Gumb_assym_log(acer_fit_X,acer_fit_Y,fin_sol(1),fin_sol(2),fin_sol(3));

%% Gumbel logistic model


[fin_sol_GL, ~, ~, sol_GL, ~] =...
    Optim_Bivar_GL(ACER_hat_mean, CI_BI, acer_fit_X, acer_fit_Y, ro(2));

Gumb_log_ACER_fit = Gumb_log(acer_fit_X,acer_fit_Y, fin_sol_GL(1)+fin_sol_GL(2));
Gumb_log_ACER_hat = Gumb_log(acer_fit_X,acer_fit_Y,1/sqrt(1 - ro(2))); % _estim

Gumb_log_MOM = Gumb_mom(bar_lev_X, bar_lev_Y, 1/sqrt(1-ro(1)), Xmom, Ymom);

% Gumb_log_MOM =  -log( Gumb_mom(bar_lev_X, bar_lev_Y, 1/sqrt(1-ro(1)), Xmom, Ymom) )/...
%     (Npoints-k_choise+1);

% % fitted assymetric logistic model:
% ReturnAL = @(x,y,fsol_X,fsol_Y,fsol) Gumb_assym_log(epsilon(x,fsol_X),epsilon(y,fsol_Y),fsol(1),fsol(2),fsol(3));
% % estimated Gumbel logistic model
% ReturnGL_estim = @(x,y,fsol_X,fsol_Y,m) Gumb_log(epsilon(x,fsol_X),epsilon(y,fsol_Y),m);
% % Gumbel logistic model by Method Of Moments
% ReturnGL_MOM = @(x,y,m,MOMx,MOMy) Gumb_mom(x,y,m, MOMx, MOMy);

%% contour fit
% close all; clear all; clc;
% load Sula_and_NORDOYAN/results_temp;

m1 = max(max(log10(ACER_hat_mean)));
m2 = max(max(log10(Gumb_assym_log_ACER_fit)));
m3 = max(max(log10(Gumb_log_ACER_hat)));
a = (min([m1, m2, m3]));
%     a = -1.8;
clear m1 m2 m3
m1 = min(min(log10(ACER_hat_mean)));
m2 = min(min(log10(Gumb_assym_log_ACER_fit)));
m3 = min(min(log10(Gumb_log_ACER_hat)));
b = (max([m1, m2, m3]))+0.1;
%     b = -3.7;
clear m1 m2 m3
v1 = linspace(a,(b+a)/2,8);
v2 = linspace((b+a)/2,b,5);
v = [v1, v2(2:end)];
clear a b v1 v2
% [X_mesh,Y_mesh] = meshgrid(bar_lev_X, bar_lev_Y);

stepX = bar_lev_X(2)-bar_lev_X(1);
blx = [bar_lev_X(1:end-1), bar_lev_X(end):stepX:bar_lev_X(ceil(numel(bar_lev_X)/2)) + bar_lev_X(end)];
clear stepX
stepY = bar_lev_Y(2)-bar_lev_Y(1);
bly = [bar_lev_Y(1:end-1), bar_lev_Y(end):stepY:bar_lev_Y(ceil(numel(bar_lev_Y)/2)) + bar_lev_Y(end)];
clear stepY

Gumb_assym_log_ACER_fit = Gumb_assym_log(epsilon(blx, fin_sol_X),epsilon(bly, fin_sol_Y),...
    fin_sol(1),fin_sol(2),fin_sol(3));
Gumb_log_ACER_fit = Gumb_log(epsilon(blx, fin_sol_X),epsilon(bly, fin_sol_Y), fin_sol_GL(1)+fin_sol_GL(2));
Gumb_log_ACER_hat = Gumb_log(epsilon(blx, fin_sol_X),epsilon(bly, fin_sol_Y),1/sqrt(1 - ro(2))); % _estim
Gumb_log_MOM = Gumb_mom(blx, bly, 1/sqrt(1-ro(1)), Xmom, Ymom);

precision = 0.05;

% v = [v(10),v(3:9),v(11:12)];
v = [v(5),v(2:4),v(6:12)];
figure;
clf
for ii = 1:numel(v) %[10,3:9,11:12]  % 3:numel(v) % 1:numel(v)   numel(v):-1:3 [6,3:5,7:12]
    
    CtableACER{ii} = contourc(bar_lev_X,bar_lev_Y, log10(ACER_hat_mean),[v(ii) v(ii)]);
    CtableACER{ii}(:,CtableACER{ii}(1,:)==v(ii)) = [];
    if ~isempty(CtableACER{ii})
        plot(CtableACER{ii}(1,:),CtableACER{ii}(2,:),'.',...
            'Color',[0 0 0],'LineWidth',2,'Markersize',10)
    end
    hold on

%     CtableACER_real{ii} = contourdiscrete2( bar_lev_X, bar_lev_Y, log10(ACER_hat_mean), v(ii), precision );
%     if ~isempty(CtableACER_real{ii})
%         plot(CtableACER_real{ii}(1,:),CtableACER_real{ii}(2,:),'.',...
%             'Color',[0 0 0]) %,'LineWidth',2,'Markersize',5)
%     end    
%     hold on
    
    CtableAL_fit{ii} = contourc(blx,bly, log10(Gumb_assym_log_ACER_fit),[v(ii) v(ii)]);
    CtableAL_fit{ii}(:,CtableAL_fit{ii}(1,:)==v(ii)) = [];
    if ~isempty(CtableAL_fit{ii})
        plot(CtableAL_fit{ii}(1,:),CtableAL_fit{ii}(2,:),'-',...
            'Color',[1 0 0],'LineWidth',1.5)
    end
    
%     CtableGL_fit{ii} = contourc(blx,bly, log10(Gumb_log_ACER_fit),[v(ii) v(ii)]);
%     CtableGL_fit{ii}(:,CtableGL_fit{ii}(1,:)==v(ii)) = [];
%     if ~isempty(CtableGL_fit{ii})
%         plot(CtableGL_fit{ii}(1,:),CtableGL_fit{ii}(2,:),'-',...
%             'Color',[0 250/255 154/255],'LineWidth',1.5)
%     end
    
%     CtableGL_hat{ii} = contourc(blx,bly, log10(Gumb_log_ACER_hat),[v(ii) v(ii)]);
%     CtableGL_hat{ii}(:,CtableGL_hat{ii}(1,:)==v(ii)) = [];
% %     if ~isempty(CtableGL_hat{ii})
% %         plot(CtableGL_hat{ii}(1,:),CtableGL_hat{ii}(2,:),'-.',...
% %             'Color',[46/255 139/255 87/255],'LineWidth',1.5)
% %     end
    
    CtableGL_MOM{ii} = contourc(blx,bly, Gumb_log_MOM,...
        exp( -10.^[v(ii) v(ii)] * (Npoints-k_choise+1) / yrs_X) );
    CtableGL_MOM{ii}(:,CtableGL_MOM{ii}(1,:)==exp( -10.^v(ii) * (Npoints-k_choise+1) / yrs_X)) = [];
    if ~isempty(CtableGL_MOM{ii}) && exp( -10.^v(ii) * (Npoints-k_choise+1) / yrs_X )>0.001
        plot(CtableGL_MOM{ii}(1,:),CtableGL_MOM{ii}(2,:),'-',...
            'Color',[65/255 105/255 225/255],'LineWidth',1.5)
    end
end

for ii = 1:numel(v)
    text(min(CtableACER{ii}(1,:)),max(CtableACER{ii}(2,:)),...
        horzcat(num2str(v(ii),'%10.2g')),...
        'HorizontalAlignment','Right','FontSize',18, 'FontName', 'Times New Roman',...
        'BackgroundColor',[0.87 0.92 0.98], 'Edgecolor',[.7 .7 .7]);
    text(max(CtableACER{ii}(1,:)),min(CtableACER{ii}(2,:)),...
        horzcat(num2str(v(ii),'%10.2g')),...
        'HorizontalAlignment','Right','FontSize',18, 'FontName', 'Times New Roman',...
        'BackgroundColor',[0.94 0.87 0.87], 'Edgecolor',[.7 .7 .7]);
end

legenda = legend( ...
    horzcat(['$\hat{\mathcal{E}}_{' num2str(k_choise) '}$']),... % horzcat(['$\hat{\mathcal{E}real}_{' num2str(k_choise) '}$']),...
    horzcat(['$\mathcal{' 'AL}_{' num2str(k_choise) '}$']),...  % horzcat(['$\mathcal{' 'GL}f' '$']),... horzcat(['$\mathcal{' 'GL}e' '$']),... 
    horzcat(['$\hat{\mathcal{' 'GL}}mm' '$']),...
    'Location',[0.525 0.9596 0.01 0.01],'Orientation','horizontal'); % [0.9105 0.893 0.001 0.002] ,'Location','NorthEastOutside'
set(legenda,'Interpreter','LaTex') % 'LaTex'
xlabel('\xi')
ylabel('\eta','Rotation',0)%'Position',[-0.095 0.5])
editplot33

clear ii precision X_mesh Y_mesh...
    CtableACER CtableACER_real CtableAL_fit CtableGL_fit CtableGL_hat CtableGL_MOM


%%
% True return level
% Return_True = @(x,y) Gumb_true(x,y,mtrue, q(1), q(2));

stepX = bar_lev_X(2)-bar_lev_X(1);
blx = [bar_lev_X(1:end-1), bar_lev_X(end):stepX:bar_lev_X(ceil(numel(bar_lev_X)/2)) + bar_lev_X(end)];
clear stepX
stepY = bar_lev_Y(2)-bar_lev_Y(1);
bly = [bar_lev_Y(1:end-1), bar_lev_Y(end):stepY:bar_lev_Y(ceil(numel(bar_lev_Y)/2)) + bar_lev_Y(end)];
clear stepY

Gumb_assym_log_ACER_fit = Gumb_assym_log(epsilon(blx, fin_sol_X),epsilon(bly, fin_sol_Y),...
    fin_sol(1),fin_sol(2),fin_sol(3));
Gumb_log_ACER_fit = Gumb_log(epsilon(blx, fin_sol_X),epsilon(bly, fin_sol_Y), fin_sol_GL(1)+fin_sol_GL(2));
Gumb_log_ACER_hat = Gumb_log(epsilon(blx, fin_sol_X),epsilon(bly, fin_sol_Y),1/sqrt(1 - ro(2))); % _estim
Gumb_log_MOM = Gumb_mom(blx, bly, 1/sqrt(1-ro(1)), Xmom, Ymom);

precision = 0.05;

figure;
clf
for ii = 1:numel(pr) % 1:numel(v)   numel(v):-1:3
    
    CtableACER{ii} = contourc(bar_lev_X,bar_lev_Y, ACER_hat_mean, [level(ii) level(ii)] );
    CtableACER{ii}(:,CtableACER{ii}(1,:) == level(ii)) = [];
%     if ~isempty(CtableACER{ii})
%         plot(CtableACER{ii}(1,:),CtableACER{ii}(2,:),'.',...
%             'Color',[0 0 0],'Markersize',10) % ,'LineWidth',2,
%     end
%     hold on

%     CtableACER_real{ii} = contourdiscrete2( bar_lev_X, bar_lev_Y, log10(ACER_hat_mean), log10(level(ii)) , precision );
%     if ~isempty(CtableACER_real{ii})
%         plot(CtableACER_real{ii}(1,:),CtableACER_real{ii}(2,:),'.',...
%             'Color',[0 0 0]) %,'LineWidth',2,'Markersize',5)
%     end    
%     hold on
    
    CtableAL_fit{ii} = contourc(blx,bly, Gumb_assym_log_ACER_fit, [level(ii) level(ii)]);
    CtableAL_fit{ii}(:,CtableAL_fit{ii}(1,:) == level(ii)) = [];
    if ~isempty(CtableAL_fit{ii})
        plot(CtableAL_fit{ii}(1,:),CtableAL_fit{ii}(2,:),'-',...
            'Color',[1 0 0],'LineWidth',1.5)
    end
    hold on
%     CtableGL_fit{ii} = contourc(blx,bly, log10(Gumb_log_ACER_fit),log10([level(ii) level(ii)]));
%     CtableGL_fit{ii}(:,CtableGL_fit{ii}(1,:)==log10(level(ii))) = [];
%     if ~isempty(CtableGL_fit)
%         plot(CtableGL_fit{ii}(1,:),CtableGL_fit{ii}(2,:),'-',...
%             'Color',[0 250/255 154/255],'LineWidth',1.5)
%     end
    
    CtableGL_hat{ii} = contourc(blx,bly, Gumb_log_ACER_hat,[level(ii) level(ii)]);
    CtableGL_hat{ii}(:,CtableGL_hat{ii}(1,:)==level(ii)) = [];
%     if ~isempty(CtableGL_hat{ii})
%         plot(CtableGL_hat{ii}(1,:),CtableGL_hat{ii}(2,:),'-.',...
%             'Color',[46/255 139/255 87/255],'LineWidth',1.5)
%     end
    
    CtableGL_MOM{ii} = contourc(blx,bly, Gumb_log_MOM,[pr(ii) pr(ii)] );
    CtableGL_MOM{ii}(:,CtableGL_MOM{ii}(1,:) == pr(ii)) = [];
    if ~isempty(CtableGL_MOM{ii})
        plot(CtableGL_MOM{ii}(1,:),CtableGL_MOM{ii}(2,:),'-',...
            'Color',[65/255 105/255 225/255],'LineWidth',1.5)
    end
    
end

for ii = 1:numel(pr)
    
    x_upper(ii) = find(abs(CtableAL_fit{ii}(2,:)-CtableGL_hat{ii}(2,:))>=0.001,1,'last');
    x_upper(ii) = CtableAL_fit{ii}(1,x_upper(ii));
    x_lower(ii) = find(abs(CtableAL_fit{ii}(1,:)-CtableGL_hat{ii}(1,:))>=0.001,1,'first');
    x_lower(ii) = CtableAL_fit{ii}(1,x_lower(ii));
    
    y_upper(ii) = find(abs(CtableAL_fit{ii}(1,:)-CtableGL_hat{ii}(1,:))>=0.001,1,'first');
    y_upper(ii) = CtableAL_fit{ii}(2,y_upper(ii));
    y_lower(ii) = find(abs(CtableAL_fit{ii}(2,:)-CtableGL_hat{ii}(2,:))>=0.001,1,'last');
    y_lower(ii) = CtableAL_fit{ii}(2,y_lower(ii));
    text(x_lower(ii), y_upper(ii)-0.1,...
        horzcat(num2str(1/(1-exp( -level(ii) * (Npoints-k_choise+1) / yrs_X)),'%10.4g')),...
        'HorizontalAlignment','Right','FontSize',18, 'FontName', 'Times New Roman',...
        'BackgroundColor',[0.87 0.92 0.98], 'Edgecolor',[.7 .7 .7]);
    text(x_upper(ii)-0.1, y_lower(ii),...
        horzcat(num2str(1/(1-exp( -level(ii) * (Npoints-k_choise+1) / yrs_X)),'%10.4g')),...
        'HorizontalAlignment','Right','FontSize',18, 'FontName', 'Times New Roman',...
        'BackgroundColor',[0.94 0.87 0.87], 'Edgecolor',[.7 .7 .7]);
end

legenda = legend( ... % horzcat(['$\hat{\mathcal{E}}_{' num2str(k_choise) '}$']),...  % horzcat(['$\hat{\mathcal{E}real}_{' num2str(k_choise) '}$']),...
    horzcat(['$\mathcal{' 'AL}_{' num2str(k_choise) '}$']),...  % horzcat(['$\mathcal{' 'GL}f' '$']),... horzcat(['$\mathcal{' 'GL}e' '$']),... 
    horzcat(['$\hat{\mathcal{' 'GL}}mm' '$']),...
    'Location',[0.525 0.9596 0.01 0.01],'Orientation','horizontal'); % [0.9105 0.893 0.001 0.002] ,'Location','NorthEastOutside'
set(legenda,'Interpreter','LaTex')
xlabel('\xi')
ylabel('\eta','Rotation',0)%'Position',[-0.095 0.5])
editplot33
clear ii precision

