clear
%%% Best Approximation Problem
%%%% generates a random feas. set  Ax=b, x>=0 data; and a random v
%%% solves min  .5norm(v-x)^2  s.t.  Ax=b, x>=0
%%% i.e., finds the projection of v onto the feas set 
%%% Various options for choosing the opt as a (nondeg) vertex
%  calls     output = semismoothsolver(A,b,v,options); for the projection
addpath(genpath('../Solvers')); 

%seed = 2;
%rngsave = rng(seed);
%   load('badexamplenov20.mat')
rngsave = rng('shuffle');
save('rngfile','rngsave');
%load('rngfile')
%ng(rngsave)
    
%%% Options
options.maxiter = 9000;
options.tol = 1e-14;
options.profile = false;
options.modprint = 1;  % number of iterations between printing in solver
options.verbose = true;
options.loadcvgfailexample = false;
options.simplex = false;
xvertex = true;   % choose v so opt. x is a vertex
polarconev = true;   % choose obj so x-v is in polar cone 
                            %(v-x in normal cone)
randxpos = 0; 
        % 0 choose random number of xpos >m; pos number is fixed number
lengthv = 100;   % multiply v to make problem harder

degen=false;
numdegen=5;

m = 300; 
n = 4000;

if options.profile
	profile clear
	profile on
end

if options.loadcvgfailexample
    %load('counterexample.mat')
    %load('cvgncefail.mat')
    load('badexamplenov20.mat')
    %v = .94*v;
    nA = norm(A,'fro');
    A = A/nA;b = b/nA;
end

if degen
    tic
    density = .8;
    rc = .01;
    xB = 10*(rand(m,1)+.1);   % nondeg BFS also ensures strict feas.
    nperm = randperm(n,m);
    nonbas = setdiff(1:n,nperm);  % x nonbasic variables so z>0
    ndegenperm = randperm(length(nperm),numdegen);
    xB(ndegenperm) = 0;       % makes nondegen BFS a degen BFS

    A(:,nonbas) = sprandn(m,n-m,density,rc);
    A(:,nperm) = sprandn(m,m,density,rc); % ensure ONE proper basis matrix
    x0 = zeros(n,1);
    x0(nperm) = xB;   % for a ***nondeg.*** BFS
    x0 = x0/norm(x0);
    nA = norm(A,'fro');
    A = sparse(A/nA);
    b = A*x0;
    fprintf('time to generate random problem %g\n',toc)
    deginds=nperm(ndegenperm);
else
    tic
    density = .005;
    rc = .01;
    xB = 10*(rand(m,1)+.1);   % nondeg BFS also ensures strict feas.
    nperm = randperm(n);
    nperm = nperm(1:m);
    nonbas = setdiff(1:n,nperm);  % x nonbasic variables so z>0
    
    %%% To guarantee a nondegenerate BFS
    A(:,nonbas) = sprandn(m,n-m,density,rc);
    A(:,nperm) = sprandn(m,m,density,rc); % ensure ONE proper basis matrix
    x0 = zeros(n,1);
    x0(nperm) = xB;   % for a ***nondeg.*** BFS
    x0 = x0/norm(x0);
    nA = norm(A,'fro');
    A = sparse(A/nA);
    b = A*x0;
    fprintf('time to generate random problem %g\n',toc)
end

fprintf('Problem has m,n,density =: %i, %i, %i \n',m,n,density)

if polarconev   % choose obj so x*-v is in polar cone
        if xvertex % ensure a vertex is opt.
		%???are we allowing for bndry of normal cone????
	       fprintf('Choosing v to ensure a NICE vertex is optimal\n')
	       options.nperm = nperm;
	       xo = x0;  % choose the constructed nice vertex
	       options.xo = xo;  %????????temp for testing
	       zo = zeros(n,1);
	       %%%%%%%to add loss of strict complementarity
	       %%%  round to get zeros
	       %temp = rand(n-m,1)+.1;
	       temp = rand(n-m,1)-.1;
	       temp = lengthv*temp.*round(temp);
	       %zo(nonbas) = lengthv*(rand(n-m,1)+.1); % original
	       zo(nonbas) = temp;
	       yo = randn(m,1);
	       v = xo -A'*yo - zo; % is in normal cone??? on bdry???
	elseif randxpos == 0 % random number of pos in zo
	       xpnum = randi(n-m,1,1); % m+xpnum els of x > 0
	       fprintf('Choosing v to ensure opt fixed x with ')
               fprintf('m+%i = %i positive xi''s \n',xpnum,m+xpnum)
	       xo = x0; 
	       randxN = randperm(length(nonbas));
	       %% add pos els to BFS xo
	       xo(nonbas(randxN(1:xpnum))) = 10*(rand(xpnum,1)+.1);
	       xo = xo/norm(xo);
	       b = A*xo;
	       xzdeg = randi(min(m,n-m)-1,1,1); % # overlapping zeros
	       zo = zeros(n,1);
	       zo(nonbas(randxN(xpnum+1:end))) = ...
	                  lengthv*(rand(length(nonbas)-xpnum,1)+.1);
	       indszpos = find(zo>0);
	       randzpos = randperm(length(indszpos));
	       xzdeg = min(20,length(randzpos));  %????
	       zo(indszpos(randzpos(1:xzdeg))) = 0; % set degeneracy
	       fprintf('xzdeg: setting %i els of zo to zero \n',xzdeg)
	       yo = randn(m,1);
	       v = xo -A'*yo - zo;  % again normal cone?? on relbdry???
	end
else
        v = randn(n,1);
        v = lengthv*v/norm(v);   % how large should v be?????
end

if options.simplex   %???????make special nonsmooth since Jk is 1by1
    tic
    d = TKPsimplex_proj(v,b);
    timesimplex = toc;

    if options.profile
        profile clear
        profile on
    end

    tic
    output = simplexsemismooth(b,v,options);
    timesemism = toc;

    if options.profile
        profile report
    end

	fprintf('times semism/simplex  %g / %g \n',timesemism,timesimplex)
	fprintf('relres semism/simplex  %g / %g \n', ...
	        norm(sum(output.xfinal)-b)/b,norm(sum(d)-b)/b)
else
    if options.profile
        profile clear
        profile on
    end
    
    tic
    output = semismoothsolver(A,b,v,options);
    timesemism = toc;
    
if options.profile
	profile report
end
	fprintf('time semism algor  %g \n',timesemism)
% 	x = output.xfinal;
% 	z = output.zfinal;
% 	y = output.yfinal;
end

