function [X_c,outs] = ...
    IRLS_LRRS(Phi,Y,d1,d2,opts)
% Implements Iteratively Reweigthed Least Squares for low-rank and
% row-sparse matrix recovery.

%   Input:  Phi = Linear operator with m rows acting on the matrix space
%                 (d1 x d2)
%            Y  = vector corresponding to measurement data
%            d1 = number of rows of matrix
%            d2 = number of columns of matrix
%          opts = struct with options, including also the model orders K1
%                 and r

%% Read options
m       = length(Y);  % Number of measurements
K       = opts.K1;    % (Maximal) Row sparsity of X
r       = opts.r;     % (Maximal) rank of X
N0      = opts.N0;  %Maximal number of iterations for the IRLS algorithm
pJ      = opts.pJ; %Non-convexity parameter pJ corresponding to the joint-sparsity inducing term
pR      = opts.pR; %Non-convexity parameter pR corresponding to the low-rank inducing term
tol     = opts.tol; %Tolerance for relative change in Frobenius for convergence (e.g., 1e-6).
mode    = opts.mode;
verbose = opts.verbose;
lambda = 1;
if isfield(opts,'saveiterates') && ~isempty(opts.saveiterates)    
    saveiterates = opts.saveiterates;
else
    saveiterates = false;
end
d_min   = min(d1,d2);
d_max   = max(d1,d2);
mu      = opts.mu;
if isfield(opts,'objectivemode') &&  ~isempty(opts.objectivemode)
    objectivemode  = opts.objectivemode;
else
    objectivemode  = 'default';
end

if isfield(opts,'tracking')
    tracking = opts.tracking;
else
    tracking = false;
end
epsmin = 1e-14;

type_mean = 'geometric';
tangent_para = 'extrinsic';

%% Preallocate some space
X           = cell(1,N0);
U_X         = cell(1,N0);
eps_1       = zeros(1,N0);
eps_2       = zeros(1,N0);
time        = zeros(1,N0);
W_1_inv_small= cell(1,N0);
sing        = cell(1,N0);
%% Iteration of the Algorithm
k=1;
first_iterate = true;
eps_2_c = Inf;
eps_1_c = Inf;
W_2_1_inv_c      = eye(d1);
W_2_2_inv_c      = eye(d2);
W_1_diag_c = ones(d1*d2,1);
W_1_inv_diag_c = ones(d1*d2,1);
W_1_small_c = ones(d1,1);
W_1_inv_small_c = ones(d1,1);
outs = struct;

speedflag = 'fast';

tic
switch mode
    case {'CG','full'}
        Phi_Phistar_op = @(Xvec) Phi_Phistar(Xvec,Phi,d1,d2);
        S_c = ones(d_max,1);
        U_X_c=W_2_1_inv_c(:,1:d1);
        V_X_c=W_2_2_inv_c(:,1:d2);
        X_vec_c = zeros(d1*d2,1);
        if strcmp(mode,'CG')
            Z_vec_c = zeros(m,1);
        end
end
singularMatrixCount = 0;
while k <= N0
    %% Calculate the solution of the minimization problem
    if k > 1
        X_c_previous = X_c;
        first_iterate = false;
    end
    if mod(k,50) == 0 
        disp(['Begin Iteration ',num2str(k),'...']);
    end
    eps_previous = eps_2_c;
    delta_previous = eps_1_c;
    lambda_previous = lambda;
    %%%%%%%%%%%%%%%%
    %% Iteration of the Algorithm

    H_cell_c = weight_op_lowrank_prec(zeros(1,d1),zeros(1,d2),S_c,1,1,...%lambda,1,...
        type_mean,0);
    H_c=H_cell_c{1};
    
    switch mode
        case 'full'
            if first_iterate
                Z_vec = (Phi*Phi')\Y;
                X_vec_c = Phi'*Z_vec;
                N_inner_c = 0;
                relres_c = 0;
            else
                if strcmp(speedflag,'fast')
                    [Hsmall,~,~,dHsmall] = weight_op_lowrank_prec(weight_op.U, ...
                        weight_op.V,weight_op.S_c,1,1,type_mean,0);
                    dHsmall{1}=dHsmall{1}(1:r_c);
                    fac = W_1_inv_eps;
                    weight_vec_c = get_weight_vec(d1,d2,Hsmall,dHsmall,tangent_para,weight_op,0);
                    Vmat = compute_Vmat_LRRS(Phi,W_1_diag_c.*W_1_inv_eps,lambda,weight_op,d1,d2,fac);
                    
                    [C] = compute_Cmat_LRRS(Phi,W_1_small_c,W_1_diag_c,lambda,weight_op,weight_vec_c,d1,d2,fac);
                    Umat = Vmat';
                    Bmatvec = (W_1_inv_eps.*(lambda./weight_op.S_c_eps+W_1_diag_c)).^(-1);
                    Amat = Phi*(Bmatvec.*Phi');
                    if C(end,end) == 0
                        X_vec_c = X_c_previous;
                        N0 = k;
                    else
                        lastwarn('', '');
                        CinvVmat = C\Vmat;
                        [warnMsg, warnId] = lastwarn();
                        if not(isempty(warnId))
                            singularMatrixCount = singularMatrixCount + 1;
%                           
                        end
                        UCinvV = Umat*CinvVmat;
                        GoodMat=Amat-UCinvV;
                        z = GoodMat\Y;
                        Phiz = Phi'*z;
                        tmp = Bmatvec.*Phiz;
    
                        CinvVmatz= C\(Vmat*z);
                     
                        tmp2 = (Bmatvec).*reshape(tangspace_to_matrixspace(CinvVmatz,weight_op),d1*d2,1);
                        X_vec_c = tmp-tmp2;
                    end
                end
                
                N_inner_c = 0;
                relres_c = 0;
            end
            X_c = reshape(X_vec_c,[d1 d2]);
    end
    %%% Calculate eps_1^2
    X_rowsums_c = sqrt(sum(X_c.^2,2)); % (d x 1) vector
    [rearr_X_c,ind_X_c] = sort(X_rowsums_c,'descend');
    X_c_allrows = X_c;
    [s_c,S_largest,ind_absE]=update_sparspara(X_rowsums_c,[],...
        [],[],rearr_X_c(K+1),d1,[],[]);
    
   switch mode
        case 'full'
            [U_X_c,singval,V_X_c] = svd(X_c);%svd(X_c,'econ');
            
        case 'CG'
            [U_X_c,singval,V_X_c] = svd(X_c);
   end

    sing_c=diag(singval);

    if strcmp(opts.epschoice,'coupled')
        eps_1_c = max(min(eps_1_c,max(rearr_X_c(K+1),sing_c(r+1)/sqrt(lambda))),epsmin);
        eps_2_c = sqrt(lambda)*eps_1_c;
    elseif strcmp(opts.epschoice,'independent')
        eps_1_prev = eps_1_c;
        eps_2_prev = eps_2_c;
        eps_1_c = max(min(eps_1_c,rearr_X_c(K+1)),epsmin);
        eps_2_c = max(min(eps_2_c,sing_c(r+1)),epsmin);
        lambda = (eps_2_c/eps_1_c)^2; % 
    elseif strcmp(opts.epschoice,'sum_coupled')
        eps_1_tilde_c_m1 = eps_1_tilde_c;
        eps_2_tilde_c_m1 = eps_2_tilde_c;
        eps_1_tilde_c = max(min(eps_1_tilde_c,rearr_X_c(K+1)),epsmin);
        eps_2_tilde_c = max(min(eps_2_tilde_c,sing_c(r+1)),epsmin);
        eps_1_c = 0.5*(eps_1_tilde_c+ eps_2_tilde_c*eps_1_tilde_c_m1/eps_2_tilde_c_m1);
        eps_2_c = eps_1_c* eps_2_tilde_c_m1/eps_1_tilde_c_m1;
        lambda = (eps_2_c/eps_1_c)^2; % 
    end
    
    [S_c,S_c_eps] = set_weight_infovec(sing_c,eps_2_c,pR,'objective_thesis');
    r_c = length(find(sing_c>eps_2_c+epsmin));
    weight_op.U = U_X_c(:,1:r_c);
    weight_op.V = V_X_c(:,1:r_c);
    weight_op.S_c = S_c(1:r_c);
    weight_op.S_c_eps = S_c_eps;
    if (k >= 2 && k <= N0 && norm(X_c-X_c_previous,'fro')/norm(X_c,'fro') < tol) || eps_1_c == epsmin ...
            || singularMatrixCount >= 1
        N0 = k;
    end
    if verbose > 1 || (verbose == 1 && (k ==1 || k == N0))
        fprintf(['%-4s k: %2d N_inner: %2d relres: %.3e delta: %.3e sing_c_r+1: %.3e eps: %.3e r_k: %1d s_k: %1d rearr_X_c: %.3e lambda: %.3e\n'], ...
                mode,k,N_inner_c,relres_c,eps_1_c,sing_c(r+1),eps_2_c,r_c,s_c,rearr_X_c(K+1),lambda)
    end
    if tracking
        switch mode
            case 'full'
                [~,singval_allrows,~] = svd(X_c_allrows,'econ');
            case 'CG'
                [~,singval_allrows,~] = svd(X_c_allrows);
        end
        sings = diag(singval_allrows);
        rankobj_before_eps(k) = eps_2_prev.^(2).*(get_rankobjective(sings,eps_2_prev,pR,'objective_thesis')...
            -log(eps_2_prev).*d_min+d_min/2);
        rankobj_after_eps(k) = eps_2_c.^(2).*(get_rankobjective(sings,eps_2_c,pR,'objective_thesis')...
            -log(eps_2_c).*d_min+d_min/2);
        sparsobj_before_eps(k)  = eps_1_prev.^(2).*(get_rankobjective(rearr_X_c,eps_1_prev,pJ,'objective_thesis')...
            -log(eps_1_prev).*d_min+d_min/2);
        sparsobj_after_eps(k)  = eps_1_c.^(2).*(get_rankobjective(rearr_X_c,eps_1_c,pJ,'objective_thesis')...
            -log(eps_1_c).*d_min+d_min/2);
       
        if mu == 0
            obj_before_eps(k) = rankobj_before_eps(k)+ ...
                sparsobj_before_eps(k);
            obj_after_eps(k) = rankobj_after_eps(k)+ ...
                sparsobj_after_eps(k);
        else
            obj_before_eps(k) = mu*(lambda_previous * rankobj_before_eps(k)+ ...
                sparsobj_before_eps(k)) + 0.5.*sum((Phi*X_c_allrows(:)-Y).^2);
            obj_after_eps(k) = mu*(lambda * rankobj_after_eps(k)+ ...
                sparsobj_after_eps(k)) + 0.5.*sum((Phi*X_c_allrows(:)-Y).^2);
        end
        sing_Hs = sing_c;

        if strcmp(mode,'CG')
            outs.relres(k) = relres_c;
            outs.N_inner(k) = N_inner_c;
            outs.resvec{k} = resvec_c;
        end
    end
    lambdas(k) = lambda;
    outs.s_largest(k) = s_c;
    outs.r_largest(k) = r_c;
    for i=d_min+1:d_max
       S_c(i)=S_c_eps;%0; 
    end
    W_1_inv_eps = eps_1_c^(2-pJ);
    for i=1:d1 
        W_1_small_c(i) = max(X_rowsums_c(i),eps_1_c)^(pJ-2);
        W_1_inv_small_c(i) = max(X_rowsums_c(i),eps_1_c)^(2-pJ);
        for j=1:d2
            W_1_inv_diag_c(i+(j-1)*d1,1) = max(X_rowsums_c(i),eps_1_c)^(2-pJ);%X_rowsums_c(i).^2+ eps_1_c.^2)^(1-pJ/2);
            W_1_diag_c(i+(j-1)*d1,1) =  max(X_rowsums_c(i),eps_1_c)^(pJ-2);
        end
    end
    if saveiterates
        X{k}=X_c;
        W_1_inv_small{k}=W_1_inv_small_c;
        sing{k} = sing_c;
    end
    eps_1(k)=eps_1_c;
    eps_2(k)=eps_2_c;
    time(k) = toc;
    k=k+1;
end
%% Tidy up the size of the relevant matrices
if tracking
    outs.rankobj_before_eps = rankobj_before_eps(1:N0);
    outs.rankobj_after_eps = rankobj_after_eps(1:N0);
    outs.sparsobj_before_eps = sparsobj_before_eps(1:N0);
    outs.sparsobj_after_eps = sparsobj_after_eps(1:N0);
    outs.obj_before_eps = obj_before_eps(1:N0);
    outs.obj_after_eps = obj_after_eps(1:N0);
end
outs.time  = time(1:N0);
if saveiterates
    outs.X      = X(1:N0);
    outs.sing    = sing(1:N0);
    outs.W_1_inv_small = W_1_inv_small(1:N0);
else
    outs.X      = X_c;
    outs.sing   = sing_c;
    outs.W_1_inv_small = W_1_inv_small_c;
end
outs.eps_1           = eps_1(1,1:N0);
outs.eps_2           = eps_2(1,1:N0);
outs.lambdas    = lambdas(1:N0);
outs.N          = N0;
outs.s_largest  = outs.s_largest(1:N0);
end

function [out] = Phi_Phistar(Xvec,Phi,d1,d2)
out=Phi*Xvec;
out=Phi'*out;
end
