function [fin_sol, fmin, pos, sol, W] = Optim_Bivar(biACER, CI, acerX, acerY, ro)

CI_plus = biACER + CI;
CI_minus = biACER - CI;

W = (log(CI_plus)-log(CI_minus)).^(-2);

pos_index = W>0;
real_index = imag(W)~=0;
if ~isempty(pos_index)
    %fprintf('Non-positive weights of the objective function changed to positive.\n')
    W = abs(W);
elseif ~isempty(real_index)
    %fprintf('Imaginary weights of the objective function changed to real.\n')
    W = abs(W);
end
% Normalizing weights to have sum(w)=1
W = W/sum(sum(W));
clear pos_index real_index condition;

%%
F_test = @(x)sum(...
        sum(...
        W.*( log(biACER) -...
        log(Gumb_assym_log(acerX,acerY,x(1),x(2),x(3))) ).^2 )); % x = [m teta fi]

sol = 100*ones(3,4);
   
if ro == 1
    ro = ro - 1e-6;
end

warning off all;
%% 1) Global 3 param LOG-level constraint optimization

F1 = F_test; 

nlnineq = @(x) [];
nlneq = @(x) [];
nonlinfcn = @(x)deal(nlnineq(x),nlneq(x));

optst_logGlobal = optimset('Display','off',...
    'Algorithm', 'interior-point', ...
    'MaxFunEvals', 10000, 'MaxIter', 10000,...
    'TolX', 1e-12, 'TolFun', 1e-12);
[sol(1,1:3), ~] = fmincon(F1,[1/sqrt(1 - ro) 0.5 0.5],[],[],[],[],...
    [-Inf 0 0],...
    [Inf 1 1],...
    nonlinfcn,optst_logGlobal);
if  any(imag(sum(sol(1,:))))
    sol(1,:) = 100*ones(1,4);
end
sol(1,4) = F_test(sol(1,1:3));
clear F1 optst_logGlobal nlnineq nlneq nonlinfcn;

%% 2) Global 3 param LOG-level Least Squares constraint optimization

F2 = @(x) sqrt(sum(W.*( log(biACER) - ...
        log(Gumb_assym_log(acerX,acerY,x(1),x(2),x(3))) ).^2 )); % x = [m teta fi]

optst_log_all = optimset('Display','off',...
    'Algorithm', 'trust-region-reflective', ...
    'MaxFunEvals', 10000, 'MaxIter', 10000,...
    'TolX', 1e-12, 'TolFun', 1e-12);
[sol(2,1:3), ~] = lsqnonlin(F2, [1/sqrt(1 - ro) 0.5 0.5], ...
    [-Inf 0 0],...
    [Inf 1 1],...
    optst_log_all);
if  any(imag(sum(sol(2,:))))
    sol(2,:) = 100*ones(1,4);
end
sol(2,4) = F_test(sol(2,1:3));
clear F2 optst_log_all;
% %% 3) Hausdorff distance 1
% 
% F3 = @(x) HausdorffDist( log(biACER),...
%     log(Gumb_assym_log(acerX,acerY,x(1),x(2),x(3))) ); 
% 
% nlnineq = @(x) [];
% nlneq = @(x) [];
% nonlinfcn = @(x)deal(nlnineq(x),nlneq(x));
% 
% optst_logGlobal = optimset('Display','off',...
%     'Algorithm', 'interior-point', ...
%     'MaxFunEvals', 10000, 'MaxIter', 10000,...
%     'TolX', 1e-12, 'TolFun', 1e-12);
% [sol(3,1:3), ~] = fmincon(F3,[1/sqrt(1 - ro) 0.5 0.5],[],[],[],[],...
%     [1 0 0],...
%     [Inf 1 1],...
%     nonlinfcn,optst_logGlobal);
% if  any(imag(sum(sol(3,:))))
%     sol(3,:) = 100*ones(1,4);
% end
% sol(3,4) = F_test(sol(3,1:3));
% clear F3 optst_logGlobal nlnineq nlneq nonlinfcn;

% %% 4) Hausdorff distance 2
% 
% F4 = @(x) HausdorffDist( zeros(size(biACER)),...
%     W.*( log(biACER) - log(Gumb_assym_log(acerX,acerY,x(1),x(2),x(3))) ).^2 ); 
% 
% nlnineq = @(x) [];
% nlneq = @(x) [];
% nonlinfcn = @(x)deal(nlnineq(x),nlneq(x));
% 
% optst_logGlobal = optimset('Display','off',...
%     'Algorithm', 'interior-point', ...
%     'MaxFunEvals', 10000, 'MaxIter', 10000,...
%     'TolX', 1e-12, 'TolFun', 1e-12);
% [sol(4,1:3), ~] = fmincon(F4,[1/sqrt(1 - ro) 0.5 0.5],[],[],[],[],...
%     [1 0 0],...
%     [Inf 1 1],...
%     nonlinfcn,optst_logGlobal);
% if  any(imag(sum(sol(4,:))))
%     sol(4,:) = 100*ones(1,4);
% end
% sol(4,4) = F_test(sol(4,1:3));
% clear F4 optst_logGlobal nlnineq nlneq nonlinfcn;
 
%%
[fmin, pos] = min(sol(:,4));
fin_sol = sol(pos,1:3);
    

end







% %% 4) Global 4 param LOG-level LS-optimization Levenberg-Marquard
% if strcmp(penalty, 'Yes')
%     F = @(x) sqrt(exp(alpha.*(x(4).^sign(x(4) - 1) - 1))) * sqrt(W).*(log(ACER)-log(x(1))+...
%         x(3)*(eta - x(2)).^x(4)); % x = [q b a c]
% else
%     F = @(x) sqrt(W).*(log(ACER)-log(x(1))+...
%         x(3)*(eta - x(2)).^x(4)); % x = [q b a c]
% end
% 
% optst_logLM_all = optimset('Display','off',...
%     'Algorithm', 'levenberg-marquardt', ...
%     'MaxFunEvals', 10000, 'MaxIter', 10000,...
%     'TolX', 1e-12, 'TolFun', 1e-12);
% [sol(3,1:4), sol(3,5)] = lsqnonlin(F, [q0 b0 a0 c0], [], [], optst_logLM_all);
% if  any(imag(sum(sol(3,:))))
%     sol(3,:) = 100*ones(1,5);
% end
% 
% clear F optst_logLM_all;
% 



