function [MGamma,gamma,VKVlist,VZKVlist,Out] = preprocessFR(Gamma,gamma,Klist,Zlist,opt)
%function [MGamma,gamma,VKVlist,VZKVlist,Out] = preprocessFR(Gamma,gamma,Klist,Zlist,opt)
%%% Preprocess data via facial reduction
%%% 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
%%%               rhoA : data rhoA used to generate the
%%%                      reduced-density-operator consrtaint (important for accuracy)
%%%               
%%% OUTPUT : MGamma,gamma : vectorized constraint system after facial reduction; 
%%%                         MGamma*Hvec(rho) = gamma
%%%          VKVlist, VZKVlist : objective data after facial reduction; see manuscript 
%%%          Out : auxiliary output
%%%                Vrho : facial vector  
%%%                MGammaorig, gammaorig,Zlist_orig,Klistorig : original data
%%% STEPS
%%% % Step 1. FR using rhoA if available; reduce size of variables
%%% % Step 2. perform facial reduction  on all of Gamma
%%% % Step 3. FR for redefining G,Z lin transf./FR for objective
%% Copyrights
% Author: Hao Hu, Haesol Im and Henry Wolkowicz 
% Email: h92hu@uwaterloo.ca, j5im@uwaterloo.ca, hwolkowi@uwaterloo.ca
if isfield(opt,'verbose')
    verbose = opt.verbose;
else
    verbose = 'true';
end
if isfield(opt,'verbose')
    tol_licols = opt.tol_licols;
else
    tol_licols = 1e-10;
end

n = length(Gamma{1});  % size of original variable rho
n2 = n^2;
mGm = length(gamma);   % number of constraints
Klist_orig  = Klist;   % cell data for objective :Hn --> Hk
Zlist_orig  = Zlist;   % cell data for objective :Hk --> Hk
MGammaorig = zeros(mGm,n2);
Gammaorig = Gamma;
gammaorig = gamma;
for ii=1:length(gamma)
    MGammaorig(ii,:) = Hvec(Gamma{ii})';
end
%%%%%%%%%%%%%find/remove trivial lin dep rows
starttriv = tic;
rankMGorig = rank(MGammaorig);
rankMGoriggorig = rank([MGammaorig gammaorig]);
if rankMGoriggorig > rankMGorig  % inconsistent original system
    fprintf('original system Gamma rho = gamma INCONSISTENT!!\n')
    keyboard
elseif mGm > rankMGoriggorig  % trivial lin dep constraints
    %%% rank with TOL = max(size(A)) * eps(norm(A))
    ranktol = max(size(MGammaorig))*eps(norm(MGammaorig));
    [~,indepid] = licols(transpose(MGammaorig),ranktol);
    MGamma = MGammaorig(indepid,:);
    gamma = gammaorig(indepid);
    Gamma = Gammaorig(indepid);
    if verbose == 2
        fprintf('removing %i trivial redundant MGamma rows;  ',...
            mGm-size(MGamma,1))
        condMGMGt = condest(MGamma*MGamma');
        fprintf('result: condest MGamma*MGamma'' = %g \n',...
            condMGMGt)
    end
    mGm = size(MGamma,1);   % new size
else
    %%% rank with TOL = max(size(A)) * eps(norm(A))
    if verbose == 2 
        fprintf('no trivial redundant MGamma rows found\n')
    end
    MGamma = MGammaorig;
end
if verbose == 2
    fprintf('End consistency; trivial lin dep check; time taken %g\n',...
        toc(starttriv))
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%end remove trivial dep rows



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 1. FR using rhoA if available; reduce size of variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
startrhoA = tic;
VrhoA = eye(n);  % for FR using rhoA
if isfield(opt,'rhoA')  % partial FR on the experimental constraints
    rhoA = opt.rhoA;
    rhoBsize = size(Gamma{1},1)/size(rhoA,1);
    [uA,dA] = eig(rhoA);
    rA = rank(dA);   % ????what tolerance????? compare to FR below???
    if rA == length(rhoA)  % NO FR
        VrhoA = eye(rhoBsize*size(rhoA,2));
    else
        VrhoA = kron(uA(:,end-rA+1:end),eye(rhoBsize));
        %%%%change Gamma data
        MGamma = zeros(length(gamma),size(VrhoA,2)^2);
        for ii = 1:length(Gamma)
            Gamma{ii} = VrhoA'* (Gamma{ii} *VrhoA);
            MGamma(ii,:) = (Hvec(Gamma{ii}))';  % real vector
            Gamma{ii} = HMat(MGamma(ii,:));   % back to Hermitian
        end
        %%% remove any (at least one) new trivial redundant constraints
        startrhoAremove = tic;
        ranktol = max(size(MGamma))*eps(norm(MGamma));
        [~,indepid] = licols(transpose(MGamma),ranktol);
        MGamma = MGamma(indepid,:);
        gamma = gamma(indepid);
        Gamma = Gamma(indepid);
        if verbose == 2 
            fprintf('removing %i red. MGamma rows after rhoA;  ',...
                mGm-size(MGamma,1))
            condMGMGt = condest(MGamma*MGamma');
            fprintf('result: condest MGamma*MGamma'' = %g;  ',...
                condMGMGt)
            fprintf('time taken = %g\n',toc(startrhoAremove))
        end
        
    end    % of rA == length(rhoA)
    nrho = size(VrhoA,2);
    if verbose == 2 
        fprintf(2,'STEP1 ')
        fprintf('FR rhoA done: dim n = %i to n_rho = %i\n',...
            n,nrho)
    end
end   % of isfield rhoA i.e. FR with rhoA
if verbose == 2
    fprintf('          end STEP 1 rhoA; time taken = %g \n',...
        toc(startrhoA))
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  END STEP 1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 2. perform facial reduction  on all of Gamma %%%%%%%
%%% after applying VrhoA
startstep2 = tic;
if verbose == 2 
    fprintf(2,'STEP2: ')
    fprintf('starting \n')
end
startFRG = tic;
[V_rho_all,Out] = FRGNPC(Gamma,MGamma,gamma,opt);

if Out.flag == -1
    fprintf('The problem is infeasible')
    return
end

if verbose == 2 
    fprintf(2,'STEP2b: ')
    fprintf(' FR for all of Gamma took %g\n',toc(startFRG))
end


if size(V_rho_all,1) > size(V_rho_all,2) %further FR is obtained
    V_rho = V_rho_all;  % final V_rho
    if verbose == 2
        fprintf(2,'STEP2c; ')
        fprintf(' vectorize; remove redund. constr. after full FR \n')
    end
    
    for ii = 1:length(Gamma)
        Gamma{ii} = V_rho'* (Gamma{ii} *V_rho);
        Gamma{ii} = (Gamma{ii}+Gamma{ii}')/2;
    end
    V_rho = VrhoA*V_rho_all;  % final V_rho
    [MGamma,gamma,Gamma] = ...
        sHvectorize_indep_map(Gamma,gamma,tol_licols);
else
    V_rho = VrhoA;
end
nrho = size(V_rho,2);



%%%% both steps for FR of Gamma are done %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verbose == 2
    fprintf(2,'STEP2 DONE: ')
    fprintf('FR; update Gamma; dim n = %i to n_rho = %i',...
        n,nrho)
    fprintf(' time taken = %g\n',toc(startstep2))
end




% Step 3. FR for redefining G,Z lin transf./FR for objective
%%%  a) get R_rho pos def. feasible for FR Gamma R_rho = gamma
%%%          ( default R_rho = I if unavailable )
%%%  b) change Klist using V_rho (i.e., create KVlist = KVlist*V_rho)
%%%  c) find V_del: range of GV(R_rho) to find V_del with equal range
%%%  d) change Klist using V_del,  smaller GUV = V_del'* GV *V_del
%%%          ( i.e., VKVlist = V_del' * KVlist *  V_del
%%%  e) compute R_delta = GUV(R_rho)
%%%  f) find V_sig: range of ZV(R_delta)
%%%  g) change Zlist using V_sig, smaller ZUV = V_sig'* ZV * V_sig
%%%          ( i.e., VZKVlist = V_sig' * Zlist*KVlist *  V_sig
%%%

if verbose == 2
    fprintf(2,'STEP3: ')
    fprintf(' start FR on obj., G,Z transf using orth \n')
end


%%% step a)  % find a pos def feas soln for Gamma(rho)=gamma conatrsint
Rhat = eye(nrho);
opt.rho0 = Rhat;   % starting point in algorithm GNQKD


%%% start on Klist
%%% step b) % first FR for operator G
KVlist = Klist;  % for Ghat operator; for VZKVlist for ZGhat operator
if size(V_rho,2) < size(V_rho,1)  % FR was done so update Klist
    for ii = 1:length(Klist)
        KVlist{ii} = Klist{ii}*V_rho;
    end
end


%%% step c) find V_del: range of GV(R_rho) to find V_del with equal range
GVofI = sum_congruence(Rhat,KVlist,0);
[ug,dg] = eig(GVofI);  %
rdg = rank(dg);

%%% step d) change Klist using V_del,  smaller GUV = V_del'* GV *V_del
VKVlist = KVlist;  % to pass to solver, i.e. this is the Ghat operator
if rdg == length(GVofI)  % no FR on f1 needed
    V_del = eye(length(GVofI));   % NOT used in this case
    if verbose == 2
        fprintf('         KVlist size %i %i not changed\n',...
            size(KVlist{1}));
    end
else
    V_del = ug(:,end-rdg+1:end);
    for ii = 1:length(KVlist)
        VKVlist{ii} = V_del'*KVlist{ii};
    end
    if verbose == 2
        fprintf('         KVlist size %i %i for f1 redefined to size %i %i \n',...
            size(Klist_orig{1}),size(VKVlist{1}));
    end
end  %  if rdg ==

%%% start on Zlist
%%% step e) compute R_delta = GUV(R_rho)
R_delta = sum_congruence(Rhat,VKVlist,0);  % after rotation of G oper

%%% step f) find V_sig: range of ZV(R_delta)
ZVofI = sum_congruence(V_del*R_delta*V_del',Zlist,0);
[uz,dz] = eig(ZVofI);
rdz = rank(dz);
V_sig = uz(:,end-rdz+1:end);

%%% step g) change Zlist using V_sig, smaller ZUV = V_sig'* ZV * V_sig
VZKVlist = cell(length(Klist_orig)*length(Zlist_orig),1);
kk = 0;
for ii = 1:length(Zlist_orig)
    for jj = 1:length(KVlist)
        kk = kk+1;
        VZKVlist{kk} = Zlist_orig{ii}*KVlist{jj};  %Z of G_V
    end
end
if rdz < length(ZVofI)  % FR on f2 needed
    for ii = 1:length(VZKVlist)
        VZKVlist{ii} = V_sig'*VZKVlist{ii};
    end
end


if verbose == 2
    fprintf('         Zlist size %i %i for f2 redefined to size %i %i \n',...
        size(Zlist_orig{1}),size(VZKVlist{1}));
end
if verbose == 2
    fprintf(2,'STEP3 ')
    fprintf(' DONE: FR for lin transf of objectives f1,f2 \n')
end

if verbose >= 1
    fprintf(2,'Facial Reduction Report \n')
    fprintf('    Reduction on gamma: %d -> %d  \n', length(gammaorig),length(gamma))        
    fprintf('    Reduction on rho  : %d -> %d  \n', length(Gammaorig{1}),length(Gamma{1}))    
    fprintf('    Reduction on delta: %d -> %d  \n', length(Zlist_orig{1}),length(VKVlist{1}))    
    fprintf('    Reduction on sigma: %d -> %d  \n', length(Zlist_orig{1}),length(VZKVlist{1}))    
end

%%
Out.V_rho = V_rho;  % facial vector 
Out.MGammaorig = MGammaorig; % origianl constraint matrix
Out.gammaorig = gammaorig;   % original rhs of the constraint system
Out.Zlist_orig = Zlist_orig; % original key map
Out.Klist_orig = Klist_orig; % original Kraus operator

end   %   of if instance  skipping FR if done already
