function [Xr,outs,opts,alg_name_out] = run_LRJS_algos(y,A,n1,n2,alg_name,opts_custom)
%run_LRJS_algos This function compares different algorithms for a given
%low-rank and joint sparse recovery problem. 
%       Inputs: y = (m x 1) measurement vector.
%               A = action of measurement operator (m x (n1*n2)) given as
%                       function handle
%              n1 = number of rows in data matrix
%              n2 = number of columns in data matrix
%        alg_name = (1 x nr_algos) cell of character strings indicating
%                    which algorithm to use.
%            opts = (1 x nr_algos) cell of custom options for the
%                   respective algorithms.
%
%               r = rank of ground truth
%              K1 = row sparsity of ground truth
%              K2 = column sparsity of ground truth
m     = length(y);

opts_extendnames = {'p','pencil_para','mode_eps','type_mean'};
alg_name_extend = 'IRLS';
[opts_new,alg_name,alg_name_out] = extend_algopts(opts_custom,opts_extendnames,...
    alg_name,alg_name_extend);

nr_algos = length(alg_name);
opts     = cell(1,nr_algos);
outs=cell(1,nr_algos);

for l=1:nr_algos
    if ~contains(alg_name{l},'IRLS')
        opts{l} = setExtraOpts(opts{l},opts_new{l});
    end
    if strcmp(alg_name{l},'RiemannianIHT_adap')
       if ~isnumeric(A)
            A = A{1}(eye(n1*n2));
       end
       if isfield(opts{l},'verbose') == false || isempty(opts{l}.verbose)
            opts{l}.verbose = true;
       end
       [Xr{l},outs{l}] = riemannian_adaptive_iht(A,n1,n2,y,opts{l});
       
   elseif contains(alg_name{l},'IRLS-LRRS')
        opts{l} = getDefaultOpts_IRLS;
        opts{l} = setExtraOpts(opts{l},opts_new{l});
        if isfield(opts{l},'K1') == false || isempty(opts{l}.K1)
            opts{l}.K1 = [];
        end
        if isfield(opts{l},'r') == false || isempty(opts{l}.r)
            opts{l}.r = [];
        end
        if isfield(opts{l},'N0') == false || isempty(opts{l}.N0)
            opts{l}.N0 = 50;
        end
        if isfield(opts{l},'pJ') == false || isempty(opts{l}.pJ)
            opts{l}.pJ = 0;     % Non-convexity parameter pJ corresponding to the joint sparsity inducing term
        end
        if isfield(opts{l},'pR') == false || isempty(opts{l}.pR)
            opts{l}.pR = 0;     % Non-convexity parameter pR corresponding to the low-rank inducing term
        end
        if isfield(opts{l},'tol') == false || isempty(opts{l}.tol)
            opts{l}.tol     = 1e-4; % %Tolerance for relative change in Frobenius for convergence (e.g., 1e-6)
        end
        if isfield(opts{l},'lambdarange') == false || isempty(opts{l}.lambdarange)
            error('Please specify a range for lambda values for the stability plot.')
        else
            lambdarange{l} = opts{l}.lambdarange;
            if length(lambdarange{l}) > 2
                 gridsize{l} = length(lambdarange{l});
            else
                if isfield(opts{l},'gridsize') == false || isempty(opts{l}.gridsize)
                    gridsize{l} = 10;
                else
                    gridsize{l} = opts{l}.gridsize;
                end
            end
        end
        if isfield(opts{l},'stabrule') == false || isempty(opts{l}.stabrule)
            opts{l}.stabrule    = 'LowRankJointSparse';
        end
        
        if ~isnumeric(A)
            A = A{1}(eye(n1*n2));
        end
        
        %        [outs{l}.x_evol,outs{l}.N,outs{l}.time,outs{l}.eps,outs{l}.diff_out,...
            %outs{l}.S_out,outs{l}.U_out,outs{l}.V_perp,outs{l}.x_evol_denoised] = ...
        if contains(alg_name{l},'IRLS-LRRS')
            lambdas = lambdarange{l};
            [Xr{l},outs_algo_c] = ...
            IRLS_LRRS(A,y,n1,n2,opts{l});
            fn = fieldnames(outs_algo_c);
            for i = 1:length(fn)
                outs{l}.(fn{i}) = outs_algo_c.(fn{i});
            end
%             if length(lambdas) == 1
%                 outs{l}.Xhat = Xr{l};
%             end
        else
            if strcmp(alg_name{l},'IRLS-maxLRRS')
            [outs{l}.X,outs{l}.eps_1,outs{l}.eps_2,outs{l}.singval,...
                outs{l}.N,outs{l}.time,...
                outs{l}.W_1_inv_diag,outs{l}.W_2_1_inv,outs{l}.W_2_2_inv,...
            outs{l}.U_X]=...
            IRLS_maxLRRS(A,y,n1,n2,opts{l});
            end
        end
        %Xr{l}=outs{l}.X{end};
        
   elseif strcmp(alg_name{l},'ATLAS')
        if isfield(opts{l},'r') == false || isempty(opts{l}.r)
            error('Please specify a rank parameter for ATLAS.')
        end
        if ~isnumeric(A)
            A = A{2}(eye(n1*n2));
        end
        if isfield(opts{l},'N0') == false || isempty(opts{l}.N0)
            opts{l}.N0 = 50;
        end
        if isfield(opts{l},'alpha') == false || isempty(opts{l}.alpha)
            opts{l}.alpha = 0.005;
        end
        if isfield(opts{l},'beta') == false || isempty(opts{l}.beta)
            opts{l}.beta = 0.005;
        end
        if isfield(opts{l},'ISTAtolerance') == false || isempty(opts{l}.ISTAtolerance)
            opts{l}.ISTAtolerance = 1e-8;
        end
        
        [u0,v0] = StartValue(y,A,n2,n1,opts{l}.r);%StartValue(y,A,n1,n2,opts{l}.r);
        
        [U,V,outs{l}.Uout,outs{l}.Vout] = ATLAS_SIMULATIONVERSION(y,A,...
            opts{l}.alpha,opts{l}.beta,u0,v0,opts{l}.N0,opts{l}.r,opts{l}.ISTAtolerance);
        Xr{l}=(U*V')';
        
        outs{l}.N=opts{l}.N0;
        for i=1:outs{l}.N
           outs{l}.X{i}= (outs{l}.Uout{i}*outs{l}.Vout{i}')';
        end
        
   elseif strcmp(alg_name{l},'SparsePowerFactorization')
        if isfield(opts{l},'r') == false || isempty(opts{l}.r)
            error('Please specify a rank parameter for SparsePowerFactorization.')
        end
        if ~isnumeric(A)
            A = A{1}(eye(n1*n2));
        end
        if isfield(opts{l},'N0') == false || isempty(opts{l}.N0)
            opts{l}.N0 = 50;
        end
        if isfield(opts{l},'alpha') == false || isempty(opts{l}.alpha)
            opts{l}.alpha = 0.5;
        end
        if isfield(opts{l},'beta') == false || isempty(opts{l}.beta)
            opts{l}.beta = 0.5;
        end
        if isfield(opts{l},'ISTAtolerance') == false || isempty(opts{l}.ISTAtolerance)
            opts{l}.ISTAtolerance = 1e-8;
        end
        if isfield(opts{l},'K1') == false || isempty(opts{l}.K1)
            error('Please specify a column sparsity parameter for SparsePowerFactorization.')
        end
        
        [u0,v0] = StartValue(y,A,n1,n2,opts{l}.r);

        [U,V,outs{l}] = rSPF(A,y,n1,n2,opts{l}.r,opts{l}.K1/n1,...
            1,u0,v0,opts{l}.N0,opts{l});
        Xr{l}=(U*V');
        
%         outs{l}.N=length(outs{l}.Uout);%opts{l}.N0;
%         for i=1:outs{l}.N
%            outs{l}.X{i}= (outs{l}.Uout{i}*outs{l}.Vout{i}');
%         end
       
   end
       
end


end

function opts = setExtraOpts(opts,opts_new)
    if ~isempty(opts_new)
        optNames = fieldnames(opts_new);
        for i = 1:length(optNames)
            optName = optNames{i};
            opts.(optName) = opts_new.(optName);
        end
    end

end

