function [ACER_hat_mean, CI_BI, bar_lev_X, bar_lev_Y] = ...
    BiACER( X, Y, eta_1_X, eta_1_Y, N_X, N_Y, k_memory, flagACER)%  k_memory, N_lev, eta_1, delta)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here


if  size(X,2)> size(X,1)
    X = X';
end
process(1,:,:) = X;

if  size(Y,2)> size(Y,1)
    Y = Y';
end
process(2,:,:) = Y;
Dim = size(process);  % [(Number of time series simulated)x(maximum number of extracted peaks)]


if size(Dim,2) == 3
    ML = zeros(Dim(3),1);
    for jj=1:Dim(3)
        ML(jj) = sum(isfinite(process(1,:,jj)));
    end
    clear jj
else
%     ML = 0; %zeros(Dim(2),1)
    ML = sum(isfinite(process(1,:)));
%     for jj=1:Dim(2)
%         ML(jj) = sum(isfinite(process(1,jj)));
%     end
%     clear jj
end


% bar_lev = linspace( min( min(min(X)) , min(min(Y)) ) - eps,...
%     max( max(max(X)) , max(max(Y)) ) + eps, N_lev);

%% Making equidistance barrier levels

bar_lev_X = linspace(min(min(X))-eps,max(max(X))+eps,N_X);%
% bar_lev_X = bar_lev;
D_eta_X = bar_lev_X(2)-bar_lev_X(1);   % Delta barrier level

bar_lev_Y = linspace(min(min(Y))-eps,max(max(Y))+eps,N_Y);
% bar_lev_Y = bar_lev;
D_eta_Y = bar_lev_Y(2)-bar_lev_Y(1);   % Delta barrier level

% clear X Y
%% Calculating ACER functions for each of the realizations
% eps_hat_X = zeros(Dim(3),N_barrier);
% eps_hat_X_ = zeros(Dim(3),N_barrier);
if size(Dim,2) == 3
    N = Dim(3);
    [A, B, G, C, C_] = deal(zeros(N_Y, N_X, N));
else
    [A, B, G, C, C_] = deal(zeros(N_Y, N_X));
    N = 1;
end
h = waitbar(0,horzcat(['ACER k =',' ',num2str(k_memory)]),'Name','Calculating ACER function...');
for zz = 1 : N
    
    N_D_eta = [(process(1,:,zz)-bar_lev_X(1))/D_eta_X + 1;...
        (process(2,:,zz)-bar_lev_Y(1))/D_eta_Y + 1];% Number of barrier levels each peak contains
    N_D_eta_floor = floor(N_D_eta);                      % Number of barrier levels each peak crosses
    N_D_eta_ceil = ceil(N_D_eta);                        % Number of barrier levels not met by each peak
    
    if k_memory==1
        for ii=1 : ML(zz)
            counter_A = 1:N_D_eta_floor(1,ii);
            counter_B = 1:N_D_eta_floor(2,ii);
            A(:,counter_A,zz) = A(:,counter_A,zz)+1;
            B(counter_B,:,zz) = B(counter_B,:,zz)+1;
            G(counter_B,counter_A,zz) = G(counter_B,counter_A,zz)+1;
        end
        %         C = ML(zz)*ones(size(A));
        C(:,:,zz) = ML(zz)*ones(N_Y, N_X);
        C_(:,:,zz) = ML(zz)*ones(N_Y, N_X); 
    else
        kk = k_memory-1;
        C_(:,:,zz) = (ML(zz)-kk)*ones(N_Y, N_X);
        for ii = kk+1 : ML(zz)
            counter_C = [max(N_D_eta_ceil(1,ii-kk:ii-1));...
                max(N_D_eta_ceil(2,ii-kk:ii-1))];
            %% NOTA BENE!!!
            counter_C(counter_C < 0) = 1;
            %%
            C(counter_C(2):N_Y, counter_C(1):N_X, zz) =...
                C(counter_C(2):N_Y, counter_C(1):N_X, zz)+1;
            
            counter_A = counter_C(1):N_D_eta_floor(1,ii);
            counter_B = counter_C(2):N_D_eta_floor(2,ii);
            
            if ~isempty(counter_A)
                A(counter_C(2):N_Y,counter_A,zz) =...
                    A(counter_C(2):N_Y,counter_A,zz)+1;
            end
            
            if ~isempty(counter_B)
                B(counter_B,counter_C(1):N_X,zz) =...
                    B(counter_B,counter_C(1):N_X,zz)+1;
            end
            
            if ~isempty(counter_A) && ~isempty(counter_B)
                G(counter_B,counter_A,zz) = G(counter_B,counter_A,zz)+1;
            end
        end
    end
    % ACER_hat(:,:,zz) = (counterA(:,:,zz)./counterC(:,:,zz) +...
    %     counterB(:,:,zz)./counterC(:,:,zz) -...
    %     counterG(:,:,zz)./counterC(:,:,zz));
    waitbar(zz/N);
end
close(h);
clear h;

C(C==0) = eps;
ACER_hat = (A./C + B./C - G./C);
ACER_hat_ = (A./C_ + B./C_ - G./C_);
clear ii jj kk zz;
clear counter_A counter_B counter_C D_eta_X D_eta_Y N_D_eta_floor N_D_eta_ceil N_D_eta
clear A C B G
%% Calculating Mean & Confidence Intervals of ACER functions (optional)
% if size(Dim,2) == 3
%     Multi_realizations = true;
if (flagACER == 1)
    ACER_hat_mean = mean(ACER_hat,3);
    ACER_hat_std = std(ACER_hat,0,3);
else
    ACER_hat_mean = mean(ACER_hat_,3);
    ACER_hat_std = std(ACER_hat_,0,3);
end

if size(Dim,2) == 3
    conf_coef = tinv((1+0.95)/2,N-1);
    CI_BI = conf_coef*ACER_hat_std/sqrt(N);     % Confidence interval
else
    conf_coef = norminv((1+0.95)/2);       % Calculating confidence interval coefficient i.e. 1.96 for 95% confidence interval
    CI_BI = conf_coef*sqrt(ACER_hat_mean)/sqrt(sum(ML)-k_memory+1);
end
% Dxy_hat_mean = diag(ACER_hat_mean)';
% Dxy_hat_std = diag(ACER_hat_std)';
% CI_Dxy = conf_coef*Dxy_hat_std/sqrt(N);
%     [fin_sol_Dxy, ~, sol_Dxy] = Opt(cat(1,squeeze(process(1,:,:)),squeeze(process(1,:,:))),...
%         bar_lev, Dxy_hat_mean, CI_Dxy, eta_1,0.9);

%% cutting
% bivariate
delta = 0.95;
[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 );

end

