function  run_plot_tests_nondeg_Semism(nn, mm, PP, numProbs, teststr, numstr, opts)
%function run_tests(r, nn, mm, NF, RR, numProbs, filename, opts)
%%% This function draws tables and plots that compare the performances of 
%%% Semism, InexactSemism, HLWB, and LSQ for varying m,n and density

addpath(genpath('..\..\Solvers'));
addpath(genpath('..\..\AuxiliaryFiles'));

%%% For Semismooth
timesemism = zeros(length(nn), numProbs);
itersemism = zeros(length(nn), numProbs);
relresidssemism = zeros(length(nn), numProbs);
reldistsolusemism = zeros(length(nn), numProbs);
mostitersemism = zeros(length(nn), 1);

%%% For the Semismooth testing solver 
timeInexactsemism = zeros(length(nn), numProbs);
iterInexactsemism = zeros(length(nn), numProbs);
relresidsInexactsemism = zeros(length(nn), numProbs);
reldistsoluInexactsemism = zeros(length(nn), numProbs);
mostiterInexactsemism = zeros(length(nn), 1);

%%% For HLWB
timeHLWB = zeros(length(nn), numProbs);
iterHLWB = zeros(length(nn), numProbs);
relresidsHLWB = zeros(length(nn), numProbs);
reldistsoluHLWB = zeros(length(nn), numProbs);
mostiterHLWB = zeros(length(nn), 1);

%%% For lsqlin
timeLSQ = zeros(length(nn), numProbs);
iterLSQ = zeros(length(nn), numProbs);
relresidsLSQ = zeros(length(nn), numProbs);
reldistsoluLSQ = zeros(length(nn), numProbs);
mostiterLSQ = zeros(length(nn), 1);

%%% Options for solvers
options.maxiter = 2000;
options.tol = opts.tol;
options.profile = false;
options.verbose = false;

for ll = 1:length(nn)
    %%% Problem parameters
    params.m    = mm(ll);
    m=params.m;
    params.n    = nn(ll);
    n=params.n;
    %params.tolerrank = max(m,n)*eps; % change below after Zpart????
    %tolerrank=params.tolerrank;
    %params.noiselevel   = NF(ll); % noise level related to the magnitude of the entries
    %noiselevel=params.noiselevel;
    params.density    = PP(ll);   % density
    p = params.density;
    params.simplex = false;
    simplex = params.simplex;
    params.xvertex = true;          % choose v so opt. x is a vertex
    xvertex = params.xvertex;   
    params.polarconev = true;       % choose obj so x-v is in polar cone
    polarconev = params.polarconev; % (v-x in normal cone)
    params.randxpos = 0; % 0 choose random number of xpos >m; pos number is fixed number
    randxpos = params.randxpos;
    params.lengthv = 0.1; % multiply v to make problem harder
    lengthv = params.lengthv;
    params.degen=false;
    degen=params.degen;
    
    %%%%% take average over numProbs (done in the end) %%%%

    maxitersemism = 0;
    maxiterInexactsemism = 0;
    maxiterHLWB = 0;
    maxiterLSQ = 0;
    for pp=1:numProbs %%% so this is to test numProbs for each test class
        
        params.seed = pp;
        rng(params.seed)
        %params.seed = 'shuffle';  % parameter for rng(seed)
       
        if opts.verbose
            disp('-------------------------------------------------');
            disp(params);
        end

	    %% generate random problem
           [A,b,v,x0,zo,options] = genProb(params,options);
           fprintf('New Problem with m,n,density =: %i, %i, %i \n',m,n,p)
    
    
        %% solve for the generated problem

    
        %%% Solve problem with Semismooth
        tic
        outputsemism = semismoothsolver(A,b,v,options);
        %output = semismooth_Sup(A,b,v,options);
        timesemism(ll,pp) = toc;
        itersemism(ll,pp) = outputsemism.iter;
        relresidssemism(ll,pp) = outputsemism.relresid(end);
        reldistsolusemism(ll,pp) = norm(outputsemism.xfinal - x0)/(1+norm(x0));
    
        if outputsemism.iter > maxitersemism
                    maxitersemism = outputsemism.iter;
        end
    
        %%% Solve problem with Semismooth testing solver
        tic
        Inexactoutputsemism = Inexactsemismoothsolver(A,b,v,1,options);
        timeInexactsemism(ll,pp) = toc;
        iterInexactsemism(ll,pp) = Inexactoutputsemism.iter;
        relresidsInexactsemism(ll,pp) = Inexactoutputsemism.relresid;
        reldistsoluInexactsemism(ll,pp) = norm(Inexactoutputsemism.xfinal - x0)/(1+norm(x0));
    
        if Inexactoutputsemism.iter > maxiterInexactsemism
            maxiterInexactsemism = Inexactoutputsemism.iter;
        end
        
        %%% Solve problem with HLWB
        [Ascaled,bscaled,D] = Diagonal_Matrix_Scaling(A,b);
        tic 
        outputHLWB = HLWBsolver(Ascaled,bscaled,v,options);
        timeHLWB(ll,pp) = toc;
        iterHLWB(ll,pp) = outputHLWB.iter;
        relresidsHLWB(ll,pp) = outputHLWB.relresid(end);
        reldistsoluHLWB(ll,pp) = norm(outputHLWB.xfinal - x0)/(1+norm(x0));
    
        if outputHLWB.iter > maxiterHLWB
            maxiterHLWB = outputHLWB.iter;
        end
    
        %%% Solve problem with lsqlin
        optionsLSQ = optimoptions(@lsqlin, 'MaxIterations',2000,'OptimalityTolerance',options.tol,...
                                            'ConstraintTolerance',options.tol, 'display', 'off');
        fprintf('\nUsing lsqlin;  NEW proj. problem on ')
        fprintf(' Ax=b,x>=0 with A %i by %i \n',m,n)
        tic
        [xfinal, ~,~,~,outputLSQ,~] = lsqlin(eye(n),v,[],[],A,b,zeros(n,1),[],max(v,0),optionsLSQ);
        timeLSQ(ll,pp) = toc;
        outputLSQ.xfinal = xfinal;
        relresid = outputLSQ.firstorderopt/(1+norm(b));
        outputLSQ.relresid = relresid;
        VarPrint = {'relres/optcond'};
        VarValue = {outputLSQ.relresid};
        fprintf('%-9s', 'iter')
        fprintf('%-20s', VarPrint{:})
        fprintf('\n')
        fprintf('%-9i', outputLSQ.iterations)
        fprintf('%-20e', VarValue{:})
        fprintf('\n')
    
        iterLSQ(ll,pp) = outputLSQ.iterations;
        try
            relresidsLSQ(ll,pp) = outputLSQ.relresid;
        catch
            relresidsLSQ(ll,pp) = 1e-2;
        end
    
        try
            reldistsoluLSQ(ll,pp) = norm(outputLSQ.xfinal - x0)/(1 + norm(x0));
        catch
            reldistsoluLSQ(ll,pp) = norm(max(v,0) - x0)/(1 + norm(x0));
        end
    
        if outputLSQ.iterations > maxiterLSQ
            maxiterLSQ = outputLSQ.iterations;
        end
    
        fprintf('Semism time: %g, InexactSemism time: %g, HLWB time: %g, LSQ time: %g \n \n', ...
                timesemism(ll,pp), timeInexactsemism(ll,pp), timeHLWB(ll,pp), timeLSQ(ll,pp))
        
    end  % end of for numProbs 

    mostitersemism(ll) = maxitersemism;
    mostiterInexactsemism(ll) = maxiterInexactsemism;
    mostiterHLWB(ll) = maxiterHLWB;
    mostiterLSQ(ll) = maxiterLSQ;
end % end of for length(nn)

%%% Organize Data for plots

for jj = 1:length(nn) 
  %%% Computing averages of Semism performance across 
    %%% i = 1:length(nn) datasets
    % semi
    semi.avertimesemi(jj) = mean(timesemism(jj,:));
    semi.averitersemi(jj) = mean(itersemism(jj,:));
    semi.averrelresidssemi(jj) = mean(relresidssemism(jj,:));
    semi.averreldistsolusemi(jj) = mean(reldistsolusemism(jj,:));

    %%% Computing medians of Semism performance across 
    %%% i = 1:length(nn) datasets
    semi.medtimesemi(jj) = median(timesemism(jj,:));
    semi.meditersemi(jj) = median(itersemism(jj,:));
    semi.medrelresidssemi(jj) = median(relresidssemism(jj,:));
    semi.medreldistsolusemi(jj) = median(reldistsolusemism(jj,:));

    semi.mostitersemi(jj) = mostitersemism(jj);
    
    %%% Computing averages of Semism testing solver performance across 
    %%% i = 1:length(nn) datasets
    % Inexactsemi
    Inexactsemi.avertimeInexactsemi(jj) = mean(timeInexactsemism(jj,:));
    Inexactsemi.averiterInexactsemi(jj) = mean(iterInexactsemism(jj,:));
    Inexactsemi.averrelresidsInexactsemi(jj) = mean(relresidsInexactsemism(jj,:));
    Inexactsemi.averreldistsoluInexactsemi(jj) = mean(reldistsoluInexactsemism(jj,:));

    %%% Computing medians of Semism testing solver performance across 
    %%% i = 1:length(nn) datasets
    Inexactsemi.medtimeInexactsemi(jj) = median(timeInexactsemism(jj,:));
    Inexactsemi.mediterInexactsemi(jj) = median(iterInexactsemism(jj,:));
    Inexactsemi.medrelresidsInexactsemi(jj) = median(relresidsInexactsemism(jj,:));
    Inexactsemi.medreldistsoluInexactsemi(jj) = median(reldistsoluInexactsemism(jj,:));

    Inexactsemi.mostiterInexactsemi(jj) = mostiterInexactsemism(jj);
    
    %%% Computing averages of HLWB testing solver performance across 
    %%% i = 1:length(nn) datasets
    % HLWB
    HLWB.avertimeHLWB(jj) = mean(timeHLWB(jj,:));
    HLWB.averiterHLWB(jj) = mean(iterHLWB(jj,:));
    HLWB.averrelresidsHLWB(jj) = mean(relresidsHLWB(jj,:));
    HLWB.averreldistsoluHLWB(jj) = mean(reldistsoluHLWB(jj,:));

    %%% Computing medians of HLWB testing solver performance across 
    %%% i = 1:length(nn) datasets
    HLWB.medtimeHLWB(jj) = median(timeHLWB(jj,:));
    HLWB.mediterHLWB(jj) = median(iterHLWB(jj,:));
    HLWB.medrelresidsHLWB(jj) = median(relresidsHLWB(jj,:));
    HLWB.medreldistsoluHLWB(jj) = median(reldistsoluHLWB(jj,:));

    HLWB.mostiterHLWB(jj) = mostiterHLWB(jj);
    
    %%% Computing averages of LSQ testing solver performance across 
    %%% i = 1:length(nn) datasets
    % LSQ
    LSQ.avertimeLSQ(jj) = mean(timeLSQ(jj,:));
    LSQ.averiterLSQ(jj) = mean(iterLSQ(jj,:));
    LSQ.averrelresidsLSQ(jj) = mean(relresidsLSQ(jj,:));
    LSQ.averreldistsoluLSQ(jj) = mean(reldistsoluLSQ(jj,:));

    %%% Computing medians of LSQ testing solver performance across 
    %%% i = 1:length(nn) datasets
    LSQ.medtimeLSQ(jj) = median(timeLSQ(jj,:));
    LSQ.mediterLSQ(jj) = median(iterLSQ(jj,:));
    LSQ.medrelresidsLSQ(jj) = median(relresidsLSQ(jj,:));
    LSQ.medreldistsoluLSQ(jj) = median(reldistsoluLSQ(jj,:));

    LSQ.mostiterLSQ(jj) = mostiterLSQ(jj);
end

%%% Save structs for plots in loadData.m

semistr = strcat('semi');
Inexactsemistr = strcat('Inexactsemi');
HLWBstr = strcat('HLWB');
LSQstr = strcat('LSQ');

%%% Save Averages
% Average CPU Time    
try 
    save(strcat(teststr,'semi',numstr,'.mat'),'-struct',semistr)
    save(strcat(teststr,'Inexactsemi',numstr,'.mat'),'-struct',Inexactsemistr)
    save(strcat(teststr,'HLWB',numstr,'.mat'),'-struct',HLWBstr)
    save(strcat(teststr,'LSQ',numstr,'.mat'),'-struct',LSQstr)
catch
    %%% You may have made a mistake in choosing your folder directory,
    %%% this will save you and let you manually try again 
    keyboard
end

end
