% run code
rand('seed',2);
randn('seed',2);

clear all
fid = fopen('nonconvex.txt','w');
repeats1 = 10; % > 1
repeats2 = 5; % > 1
tol = 1e-5;
r = 1e8;

for m = [3000 5000]
    % set parameters
    n = m/5;
    k = n/5;
    
    for p = [0.5 0.6 0.7]  % p*m is the size of feasible system
        tols = m*tol; % a heuristic
        
        Vdisplay = [];
        for ii = 1:repeats1
            
            A1 = randn(m,n);% m>>n
            s0 = min(max(randn(k,1),-r),r);
            s = zeros(n,1);
            I = randperm(n);
            s(I(1:k)) = s0; % generate a sparse solution
            btmp = A1*s;
            b1 = zeros(m,1);
            b1(1:p*m) = btmp(1:p*m) + 1e-2*rand(p*m,1);
            b1(p*m+1:m) = btmp(p*m+1:m) - 50*rand(m-m*p,1); % heuristics attempting to destroy feasibility
            
            A2 = randn(m,n);
            b2 = A2*s - 50*rand(m,1); % heuristics attempting to destroy feasibility
            
            
            Vtmp = [];
            for iii = 1: repeats2
                % generate random points
                x0 = randn(n,1); % random initial point
                
                tmp = x0; %%%%%%%%%%%%%% project to be feasible
                tmp0 = max(min(tmp,r),-r);
                tmp2 = tmp.^2 - (tmp0 - tmp).^2;
                [~,I] = sort(tmp2,'descend');
                x0 = zeros(n,1);
                x0(I(1:k)) = tmp0(I(1:k));
                xtest = x0;
                
                % log penalty
                type = 'log';
                epsilon = 0.9;
                tol_init = tol;
                iter = 0; time = 0;
                % profile on
                while epsilon > tol/10
                    fprintf(' Log: inner loop: epsilon = %g\n', epsilon)
                    [x,time1,iter1] = proj_ls_nonconvex(A1,A2,b1,b2,r,k,epsilon,type,tol_init,x0,inf,0,500);
                    epsilon = epsilon/10;
                    tol_init = max(tol_init/3,tol/100);
                    x0 = x;
                    time = time + time1; iter = iter + iter1;
                end
                
                % profile off
                % compute feasibility
                xlog = x;
                
                %fprintf(' Termination: log penalty: epsilon %g iter %10d time %5.1f feas %2.1e\n\n', epsilon*10, iter, time, feas1)
                
                % errors for the feasible part
                v1 = min([A1*xlog - b1, A2*xlog - b2],[],2);
                v3 = min([A1*xtest - b1, A2*xtest - b2],[],2);
                
                size_log = sum(max(v1,0) < tols/4);
                size_log0 = sum(max(v1,0) <= 0);
                size_ran = sum(max(v3,0) < tols/4);
                size_ran0 = sum(max(v3,0) <= 0);
                
                % Save data for one problem
                Vtmp = [Vtmp; iter time size_ran0 size_ran size_log0 size_log];
            end
            % take max
            Vdisplay = [Vdisplay; mean(Vtmp(:,1:2)) max(Vtmp(:,3:6))/m];
        end
        % average and print
        fprintf(fid,'%6.0d  & %2.1f & %6.0f & %7.1f & %3.2f & %3.2f & %3.2f & %3.2f\n', m, p, mean(Vdisplay));
    end
end
fclose(fid);