function [ubd,lbd,result] = GaussNewton(Gamma,gamma,Klist,Zlist,opt)
%function [ubd,lbd,Out] = GaussNewton(Gamma,gamma,Klist,Zlist,opt)
%%% Run the Gauss-Newton approach for QKD problem
%%% step 1 : preprocess the data
%%% step 2 : run GNQKD (gauss-newtwon algorithm)
%%% INPUT : Gamma, gamma : constraint data <Gamma_i,rho>=gamma_i
%%%         Klist, Zlist : obejective data; Kraus operator, Key map
%%%         opt : struct, option
%%%               verbose : print auxiliary progress report
%%%                         0=print nothing
%%%                         1=print some outputs 
%%%                         2=print all outputs 
%%%               rhoA : data rhoA used to generate the
%%%                      reduced-density-operator consrtaint (important for accuracy)
%%%               tolerFR : facial reduction tolerance
%%%               tolerGN : gauss-newton tolerance 
%%%               tol_licols : licols tolerance 
%%% OUTPUT : ubd,lbd : best upper bound and best lower bound
%%%          Out : auxiliary output
%%%                rho: primal variable of the reduced model
%%%                Z: dual variable w.r.t. psd constraint
%%%                y: dual variable w.r.t. affine constraint
%%%                ubs: history of upper bounds
%%%                lbs: history of lower bounds
%%%                iter: total number of interation from GNQKD
%%%                V_rho: facial vector from facial reduction
%%%                origobjval: objective value 
%%%                VrhoV: V_rho*rho*V_rho'; feasible point to the original problem
%%%                gap: ubd-lbd; nomimal gap between lower and upper bounds
%%%                relgap: relative gap between lower and upper bounds
%%%                time: running time 
if isfield(opt,'verbose')
    verbose = opt.verbose;
else
    verbose = 1;
end
if isfield(opt,'rhoA')
    rhoA = opt.rhoA;  % matrix used to generate the reduced-density-operator constraints
end

start_time = tic;

Klist_orig = Klist;
Zlist_orig = Zlist;

if verbose == 2
    data_check(Gamma,gamma,Klist,Zlist)  % checking data eligibility
end

if verbose >=1
    fprintf('\nInitial sizes: Gamma{1}: n=%i , Klist{1}: %i-by-%i, Zlist{1}: %i-by-%i \n', ...
        length(Gamma{1}),size(Klist{1}),size(Zlist{1}))
end
%% STEP 1 preprocessing (facial reduction)
if verbose >= 1
    fprintf('Starting preprocessing for facial reduction\n')
end
[MGamma,gamma,VKVlist,VZKVlist,OutFR] = preprocessFR(Gamma,gamma,Klist,Zlist,opt);
MGammaorig = OutFR.MGammaorig;
gammaorig = OutFR.gammaorig;

%% Check for trivial cases; unique solution, infeasible problem
if size(MGamma,1) >= size(MGamma,2)  % m>=n^2
    if verbose >=1
        fprintf('\n In run_GNQkD: Trivial case unique solution rho! \n')
    end
    rho = HMat(MGamma\gamma);
    if min(eig(rho)) >= 0
        result.rho = rho;
        result.V_rho = OutFR.V_rho;
        VrhoV = OutFR.V_rho*rho*OutFR.V_rho';
        result.VrhoV = VrhoV;
        %calculate the current objective
        MDelta = sum_congruence(rho,VKVlist,0);
        MSigma = sum_congruence(rho,VZKVlist,0);
        EigDelta = eig(MDelta);
        EigSigma = eig(MSigma);
        objval = real(sum(EigDelta.*log(EigDelta)) - ...
            sum(EigSigma.*log(EigSigma)));
        result.objval = objval;
        result.bestub = objval;
        result.bestlb = objval;
        result.time = 0;
        result.iter = 0;

        ubd = objval;
        lbd = objval;
        if verbose >=1
            fprintf('               objective value is %g\n',...
                objval )
        end
    else
        fprintf('infeasible problem rho NOT psd  \n')
    end  % end of  min(eig(rho)) >= 0
    return
end % of if ......  % m>=n^2

%% STEP 2 call the Gauss-Newton
[rho,Out] = GNQKD(MGamma,gamma,VKVlist,VZKVlist,opt);

% save the results
result.rho = rho;
result.Z = Out.Z;
result.y = Out.y;
result.ubs = Out.ubs; % history of upper bounds
result.lbs = Out.lbs; % history of lower bounds 
result.iter = Out.noiter;
result.V_rho = OutFR.V_rho; % facial vector

% obtain the objective value by recovering V_rho*rho*V_rho' from rho
V_rho = OutFR.V_rho;
VrhoV = V_rho*rho*V_rho';
GVrhoV = sum_congruence(VrhoV,Klist_orig,0); % G(rho)
ZGVrhoV = sum_congruence(GVrhoV,Zlist_orig,0); % Z(G(rho))
GVrhoV = (GVrhoV+GVrhoV')/2; % symmetrize
ZGVrhoV = (ZGVrhoV+ZGVrhoV')/2;
% remove 0 eigenvalues to compute logm without warning
[Udel,Ddel] = eig(GVrhoV,'vector');
inddel = find(Ddel>0);
  %GVrhoV = (Udel(:,inddel).*Ddel(inddel)')*Udel(:,inddel)'; % G(rho)
f1 = Ddel(inddel)'*log(Ddel(inddel)); % compute trace(G(rho)logG(rho))
[Usig,Dsig] = eig(ZGVrhoV,'vector');
indsig = find(Dsig>0);
  %ZGVrhoV = (Usig(:,indsig).*Ddel(indsig)')*Udel(:,indsig)';
Lams = Ddel(inddel)*log(Dsig(indsig))'; % eigenvalue coefficients
Utemp = Udel(:,inddel)'*Usig(:,indsig);
Utemp = real(Utemp).^2 + imag(Utemp).^2 ;
f2 = sum(sum(Lams.*Utemp)); % compute trace(G(rho)logZ(G(rho)))

result.origobjval = f1-f2;  % objetive value 
result.VrhoV = VrhoV;  % feasible point to the orignal problem V_rho*rho*V_rho'


result.gap = result.origobjval - Out.bestlb; % gap between upper and lower bound
result.relgap = 2*(result.gap)/ ...
    (result.origobjval + Out.bestlb + 1);  % relative gap

result.time = toc(start_time);
if verbose == 2
    fprintf('\n')
    fprintf(2,'Output from GaussNewton:  \n')
    fprintf('solver time %g \n',result.time)
    relprimalfeas = norm(MGammaorig*Hvec(VrhoV)-gammaorig)/norm(gammaorig);
    fprintf('relative residual ||MGamma*VrhoV-gamma||/||gamma|| = %g \n', relprimalfeas)
    fprintf('relpsd: mineig(rho)/norm(rho) = %e, mineig(Z)/norm(Z) = %e \n',...
           min(eig(rho))/norm(rho,'fro'), min(eig(Out.Z))/norm(Out.Z,'fro'))
end

ubd = result.origobjval;    % upper bound 
lbd = Out.bestlb;           % lower bound
if verbose >= 1
    fprintf('lower bound = %g, upper bound = %g\n',lbd,ubd)
end
end
