function  run_plot_tests_LP(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 linprog for varying m,n and density

addpath(genpath('..\..\'));
addpath(genpath('C:\Program Files\MATLAB\R2022b\toolbox\optim\optim'))

%%% For Semismooth
timesemism = zeros(length(nn), numProbs);
relresidssemism = zeros(length(nn), numProbs);

%%% For linprog dual simplex
timelinprogds = zeros(length(nn), numProbs);
relresidslinprog = zeros(length(nn), numProbs);

%%% For linprog interior point
timelinprogds = zeros(length(nn), numProbs);
relresidslinprog = zeros(length(nn), numProbs);

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

for ll = 1:length(nn)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%set options
    params.m    = mm(ll);
    m=params.m;
    params.n    = nn(ll);
    n=params.n;
    params.p    = PP(ll);   % density
    p = params.p;
    maxiter = 1200;
    dexp = 14;
    nondegvertex=true;   % true for now as other generation not working yet
    tolorig = 10^(-dexp);
    setprofile = false;

    if setprofile
        profile clear
        profile on
    end 

    modprint = 1;
    lengthv=0.1;    % difficulty of problem
    %%%%%%%%%
    options.nondegvertex = nondegvertex;
    options.dexp = dexp;
    options.tol = tolorig;
    options.maxiter = 1200;  % for solver
    options.tolLP = tolorig;
    options.profile = setprofile;
    options.modprint = modprint;   % number of fprints in while in solver
    options.lengthv = lengthv;   % larger makes problem harder
    
    maxitersemism = 0;
    maxiterlinprog = 0;
    for pp=1:numProbs %%% so this is to test numProbs for each test class
    
        params.seed = pp;
        rng(params.seed)
        
        if opts.verbose
            disp('-------------------------------------------------');
            disp(params);
        end
    
        %% generate random problem
        %%%% generate problem with a nondegerate vertex as opt. (generically)
        if nondegvertex  % a nondegenerate vertex generically
            startgen = tic;
            rc = .01;  % for sparsity recip. cond. #
            xB = 10*(rand(m,1)+.1);   % nondeg BFS also ensures strict feas.
            nperm = randperm(n);
            mperm = nperm(1:m);   % indices of basic columns
            nonbas = setdiff(1:n,mperm);  % x nonbasic variables so z>0
            try
                A = spalloc(m,n,(m*p)*n);
            catch
                keyboard
            end
            A(:,nonbas) = sprandn(m,n-m,p,rc);
            %% the following, generically, makes basis matrix nonsingular
            A(:,mperm) = sprandn(m,m,p,rc); % ensure ONE proper basis matrix
            x0 = zeros(n,1);
            x0(mperm) = xB;   % for a ***nondeg.*** BFS
            x0 = x0/norm(x0);  % optimal xo norm = 1
            nA = norm(A,'fro');
            A = sparse(A/nA);
            b = A*x0;
            xo = x0;  % choose the constructed nice vertex
            zo = zeros(n,1);
            zo(nonbas) = lengthv*(rand(n-m,1)+.1);
            yo = randn(m,1);
            c = A'*yo - zo;  % dual str. feas.
            trueoptval = c'*xo;  % true optimal value
            endgentime = toc(startgen);
            fprintf('time to generate random problem %g\n',endgentime)
        else  % not polished yet??????????????
            A = randn(m,n);
            A = A/norm(A,'fro');
            x0 = [(rand(m,1)+.1) ;zeros(n-m,1)];
            x0 = x0/norm(x0);
            z0 = [zeros(m,1) ;rand(n-m,1)+.1];
            z0 = 2*z0/norm(z0);
            inds = randperm(n);
            x0 = x0(inds);
            z0 = z0(inds);
            b = A*x0;
            y0 = randn(m,1);
            y0 = 2*y0/norm(y0);
            c = A'*y0-z0;         %      A'y >= c;  c = A'y-z
        end


        %% solve for the generated problem

        %%% Solve LP with Exact Semismooth
        
        Roptvals = [];
        Roptvalsdiff = [];
        options.yk = zeros(m,1);
        
        if options.verbose
              VarPrintint = {'m','#xj>0'};
              VarPrint = {'optvalue','relresid'};
              fprintf('\nUsing nonsmooth Newton  to solve a NEW LP problem \n')
        end
        %relerroropt = inf; %%% This needs to get reassigned in while loop???
                           %%% otherwise iter grows to 500 every time
        %Rstep = 1;
        Rchange = .1;   % small increase to force new basis for warm start
        %R = 10;   % R start??? more adaptive???
        R = min(50,sqrt(m*n)*norm(b)/(norm(c)+1)); %??why use this estimate??
        Rs = R;
        iter = 0;
        countbasic = [];
        countxfinalpos = [];
        %vR = R*c;
        %vs = vR;
        bR = b/R;
        lbnd = -inf;  % lower bound on optimal value
        newlbnd = lbnd;
        ubnd = inf;
        newubnd = inf;
        cxcnt = 0;
        cycnt = 0;
        relresidprojs = inf;
        
        %%% How does R impact performance??
        tols = [];
        minoriters = [];
        normc = norm(c)+1;
        relpfdf = inf;
        startticProj = tic;   % start time for LP loop
        while (R < inf) 
            iter = iter + 1;  % major iterations
            %%% increase accuracy at each iteration for projection
            options.tol = max(tolorig, ...
            min([1e-12 .01*relpfdf .01*relresidprojs(end)])); % increase accuracy 
            tols = [tols options.tol];
	        %options.Rs = Rs;   % why??? needed??
            if options.verbose
                fprintf('\n\n LP solve starting major iter. %i \n',iter)
	         %   norm(x)/R %g ,norm(z)/R %g ,' ...
	         %         iter,norm(x)/R,norm(z)/R)
	         %  fprintf('(norm(A''*(-y/R)-c-(z)/R))/norm(c)  %g %g %g\n',...
                 %  (norm(A'*(-y/R)-c-(z)/R))/norm(c))
            end
            %output = semismoothsolver(A,b,vR,options);
            output = semismoothsolver(A,bR,c,options);
            relresidprojs = [relresidprojs output.relresid(end)];
            minoriters = [minoriters output.iter];
            xapprox = R*output.xfinal;
            %yapprox = output.yfinal;
            %zapprox = output.zfinal;
            %pfapprox = norm(A*xapprox-b);
            %dfapprox = norm(c+A'*(-yapprox)+zapprox);
            %csapprox = sum(prod([xapprox zapprox],2));
            %KKTapprox = pfapprox+dfapprox+csapprox;
            newlbnd = c'*xapprox;  %  lower bound
            %% change to no improv in lower or upper bounds?????
            if newlbnd <= lbnd   % no improvement in lower feas bnd
	        cxcnt = cxcnt+1;
            else
                lbnd = newlbnd;
	            lbnds(iter) = lbnd;
	            cxcnt = 0;
            end
        
            Roptvals = [Roptvals newlbnd];
            %relerroropt = abs(trueoptval - Roptvals(end))/(1+abs(trueoptval));
        
            %    if options.verboseplot
            %        %%% Plot optimal values
            %	    figure(1)
            %	    clf
            %	    plot([Roptvals trueoptval],'-.')
            %	    title('improving lower bounds to cvxoptval')
            %    
            %        %%% Plot difference to cvxoptval
            %        %Roptvalsdiff = [Roptvalsdiff relerroropt];
            %	%    figure(2)
            %	%    clf
            %	%    semilogy(Roptvalsdiff,'-.')
            %	%    title('relative difference between cvxoptval and Roptval')
            %    end
            
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%do Tapia indices work without quadr cnvgnce???????????
            %%%    %%% Tapia Indices : xtapia = x^{k+1} / x^k
            %%%    if iter >= 2
            %%%
            %%%        xtapia = output.xfinal./xpre;
            %%%
            %%%        %%% To handle dividing by zero
            %%%            xprezero = find(xpre == 0);
            %%%            xcurrentzero = find(output.xfinal == 0);
            %%%            
            %%%            %%% To handle NaN values:
            %%%            % We use the intersection of indices as 
            %%%            % find(xpre == 0) is not always a subset of find(output.xfinal == 0)
            %%%            xtapia(intersect(xprezero, xcurrentzero)) = 0;
            %%%        %%% To handle dividing by zero
            %%%
            %%%        %%% Based on change in x, we have the following pos. variables
            %%%        % measured by the rate of change between iterations
            %%%        xbasic = find(xtapia > 0.1);
            %%%        countbasic = [countbasic length(xbasic)];
            %%%    end
            %%%    xpre = output.xfinal;
            %%%    %%% end of Tapia Indices
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            xfinalpos = sum(output.xfinal>0);
            countxfinalpos = [countxfinalpos xfinalpos];
        
            if options.verbose
                fprintf('%-15s', 'R')
                fprintf('%-9s', VarPrintint{:});
                fprintf('%-15s', VarPrint{:});
                fprintf('\n')
                fprintf('%-15e',R) 
                fprintf('%-9i', m,xfinalpos)
                fprintf('%-15e',Roptvals(end),output.relresid(end))
                fprintf('\n')
            end
            %x = output.xfinal;
            %y = output.yfinal;
            %z = output.zfinal;
        
            %if options.verboseplot
        	%    xplot % plot sorted positive x values
            % end
            
        
            %%%% calculate new y for warm start
            wopt = output.xfinal;
            zopt = output.zfinal;
            yopt = output.yfinal;
            indsN = find(zopt > 0);
            indsB = setdiff(1:n,indsN);
            wBopt = wopt(indsB);
            zNopt = zopt(indsN);
            AB = A(:,indsB);  % basic type cols
            AN = A(:,indsN);  % nonbasic type cols
            deltaypre = lsqminnorm((AB*AB'),b);
            %bB = AB'*lsqminnorm((AB*AB'),b);
            bB = AB'*deltaypre;
            etemp = bB-R*wBopt;
            ftemp = R*bB;
            indsratio = find(etemp > 1e-6);
            if isempty(indsratio)
	            wBstep = inf;
            else
                    wBstep = min(ftemp(indsratio)./etemp(indsratio));
            end

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %bN = AN'*lsqminnorm((AB*AB'),b);
            bN = AN'*deltaypre;
            etemp = -(bN+R*zNopt);
            ftemp = -R*bN;
            indsratio = find(etemp > 1e-6);

            if isempty(indsratio)
	            zNstep = inf;
            else
                zNstep = min(ftemp(indsratio)./etemp(indsratio));
            end

            Rnew = min(wBstep,zNstep);
            if Rnew < R
	            fprintf('Rnew < R???error????\n')
	            keyboard
            elseif Rnew == inf
	            fprintf('\nopt. LP basis found in while loop, Rnew=inf; ')
                    yoptLP = lsqminnorm((A(:,indsB)'),c(indsB));
	            xoptLP = zeros(n,1);
                    xoptLP(indsB) = lsqminnorm(A(:,indsB),b);
	            %yoptLP = (A(:,indsB)')\c(indsB);
	            zoptLP = zeros(n,1);
	            zoptLP(indsN) = A(:,indsN)'*yoptLP - c(indsN);
	            fprintf(' and min zoptLP dual opt check is %g \n',min(zoptLP))
	            R = Rnew;
	           % keyboard  % ???? works?? returns where??? fix logic???
            else
	            R = Rnew;
            end
        
            if R < inf   % opt not found yet
                %options.yk = (R/(R + Rstep))*output.yfinal;   % for input for next run
                %options.yk = 0*options.yk;
                %R = R + Rstep;
                deltay = ( (1/R) - (1/Rs(end)))*deltaypre;
                deltaw = AB'*deltay;  % confirms dual feas. by construction
                deltaz = -AN'*deltay;
                pfcheck = norm(AB*deltaw - ( (1/R) - (1/Rs(end)))*b);
                pfdfcheck = -min([0;wBopt+deltaw;zNopt+deltaz]);
                if max(pfdfcheck,pfcheck) > 1e-8
	                fprintf('Rn step too large???\n')
	                keyboard
                end
                options.yk = yopt + deltay;
                R = R + Rchange/(iter^2);  % increase R slightly beyond last opt.
                Rs = [Rs R];
                bR = b/R;
            end  % of if R == inf
        end   % end of while for Rs

        stopticProj = toc(startticProj);  % end time for LP loop 
        %output = semismooth_Sup(A,b,v,options);

        pfproj = norm(A*xoptLP-b)/(1+norm(b));
        dfproj = norm(zoptLP - A'*yoptLP+c)/(1+norm(c));
        
        timesemism(ll,pp) = stopticProj;
        relresidssemism(ll,pp) =  pfproj + dfproj + xoptLP'*zoptLP;

        %%% Solve problem with linprog dual simplex
        fprintf('\nUsing linprog with dual simplex;  NEW proj. problem on \n')

        %%default algorithm is dual simplex
        linprogoptionsds = optimoptions('linprog','Display','none', 'OptimalityTolerance',1e-10);
        

        startticlinprogds = tic;
        [xlinprogds,flinprogds,flaglinprogds,outputlinprogds,lambdalinprogds] = ...
                linprog(-c,[],[],A,b,zeros(n,1),[], linprogoptionsds);
        stopticlinprogds = toc(startticlinprogds);  % end time for linprog
        
        try
            pflinprogds(ll,pp) = norm(A*xlinprogds-b)/(1+norm(b));
            dflinprogds(ll,pp) = norm(lambdalinprogds.lower - A'*lambdalinprogds.eqlin+c)/(1+norm(c));
            cslinprogds(ll,pp) = norm(lambdalinprogds.lower'*xlinprogds)/ ...
                                    (1 + max(norm(lambdalinprogds.lower),norm(xlinprogds)));
    
            timelinprogds(ll,pp) = stopticlinprogds;
            relresidsds(ll,pp) = pflinprogds(ll,pp) + dflinprogds(ll,pp) + cslinprogds(ll,pp);
        catch 
            if isempty(xlinprogds)
                pflinprogds(ll,pp) = Inf;
                dflinprogds(ll,pp) = Inf;
                cslinprogds(ll,pp) = Inf;
                timelinprogds(ll,pp) = Inf;
                relresidsds(ll,pp) = Inf;
            else
                keyboard
            end
        end

        

        %%% Solve problem with linprog interior point 
        linprogoptionsipm = optimoptions('linprog','Display','none', ...
          	'Algorithm','interior-point','OptimalityTolerance',1e-10);

        fprintf('\nUsing linprog with interior point method;  NEW proj. problem on \n')
        startticlinprogipm = tic;
        [xlinprogipm,flinprogipm,flaglinprogipm,outputlinprogipm,lambdalinprogipm] = ...
	               linprog(-c,[],[],A,b,zeros(n,1),[],linprogoptionsipm);
        stopticlinprogipm = toc(startticlinprogipm);  % end time for linprog
        
        try
            pflinprogipm(ll,pp) = norm(A*xlinprogipm-b)/(1+norm(b));
            dflinprogipm(ll,pp) = norm(lambdalinprogipm.lower - A'*lambdalinprogipm.eqlin+c)/(1+norm(c));
            cslinprogipm(ll,pp) = norm(lambdalinprogipm.lower'*xlinprogipm)/...
                (1 + max(norm(lambdalinprogipm.lower),norm(xlinprogipm)));

            timelinprogipm(ll,pp) = stopticlinprogipm;
            relresidsipm(ll,pp) = pflinprogipm(ll,pp) + dflinprogipm(ll,pp) + cslinprogipm(ll,pp);
        catch 
            if isempty(xlinprogipm)
                pflinprogipm(ll,pp) = Inf;
                dflinprogipm(ll,pp) = Inf;
                cslinprogipm(ll,pp) = Inf;
                timelinprogipm(ll,pp) = Inf;
                relresidsipm(ll,pp) = Inf;
            else
                keyboard
            end
        end
        
        
        fprintf('\nSemism time: %g, linprog dual simplex time: %g, linprog interior point time: %g \n \n', ...
                timesemism(ll,pp), timelinprogds(ll,pp), timelinprogipm(ll,pp))
            
    end % of for pp = 1:numProbs

end % of for ll = 1:length(nn)

%%% Save structs for plots in loadData.m
semi.timesemi = timesemism;
semi.relresidssemi = relresidssemism;

linprogds.timelinprogds = timelinprogds;
linprogds.relresidslinprogds = relresidsds;
linprogds.pflinprogds = pflinprogds;
linprogds.dflinprogds = dflinprogds;
linprogds.cslinprogds = cslinprogds;

linprogipm.timelinprogipm = timelinprogipm;
linprogipm.relresidslinprogipm = relresidsipm;
linprogipm.pflinprogipm = pflinprogipm;
linprogipm.dflinprogipm = dflinprogipm;
linprogipm.cslinprogipm = cslinprogipm;

%%% Save structs for plots 

semistr = strcat('semi');
linprogdsstr = strcat('linprogds');
linprogipmstr = strcat('linprogipm');
%%% Save Averages
% Average CPU Time   

try 
    save(strcat(teststr,'semi',numstr,'.mat'),'-struct',semistr)
    save(strcat(teststr,'linprogds',numstr,'.mat'),'-struct',linprogdsstr)
    save(strcat(teststr,'linprogipm',numstr,'.mat'),'-struct',linprogipmstr)
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
