%  as of/starting Nov30/15 ... polishing dec5/15
%  solves  x'Ax+2a'x  s.t. x'x<=delta   b'x<=beta
%  This code calls Algoritm3 to solve the second class of test problems and make
%  a comparsion with the SOCP/SDO reformulation
% 
%%======================================================================
%  Initialization below:
%%  notests
%%  n  for dimension
%%  density  for density of A
clear all         
clc
profile clear
profile on
format long e
global A a b beta n    
%===================================================================
%n=input('Enter n:');
%density=input('Enter density:');
 rng('shuffle')
 rngsave=rng;
 rng=rngsave;
density=.1;
ns=100:100:400;

TIMEs=zeros(1,length(ns));
TIMESTRONGs=zeros(1,length(ns));
KKT1s=zeros(1,length(ns));
KKT2s=zeros(1,length(ns));
ACCURACYNEWs=zeros(1,length(ns));
ACCURACYSOCPs=zeros(1,length(ns));
TIMESOCPs=zeros(1,length(ns));
nLNGMs=zeros(1,length(ns));
for nn=1:length(ns);
n=ns(nn);

%===================================================================
notests=10;   
KKT1=zeros(notests,1);   
KKT2=zeros(notests,1);   
STRONG=zeros(notests,1);
nLNGM=zeros(notests,1);
TIME=zeros(notests,1);
TIMESTRONG=zeros(notests,1);
TIMESOCP=zeros(notests,1);
ACCURACYNEW=zeros(notests,1);
ACCURACYSOCP=zeros(notests,1);

for ii=1:notests,
    
    fprintf('\n\n Start new problem #%i in secondclass.m with n=%i density=%g \n',ii,n,density)
    %% generate test problem 
    [A,a,b,beta,delta]=secondclass(n,density);

    fprintf('solve min x''Ax +2a''x s.t. ||x||^2<=delta=%g, b''x<=beta = %g \n',delta,beta)
    [kkt1,kkt2,lambda_1opt,strong,LNGM,time,timestrong,fval,xopt,xopt1]=Algorithm3(A,a,b,beta,delta,n);
    KKT1(ii)=kkt1;
    KKT2(ii)=kkt2;
    STRONG(ii)=strong; % if strong=0, we dont have strong duality. if strong=1, we have strong duality for eTRS
    nLNGM(ii)=LNGM;  % LNGM=0, solution of eTRS does not happen at LNGM of TRS. if LNGM=1,solution of eTRS is LNGM of TRS.         
    TIME(ii)=time;   % Time taken for Algorithm 3 (Main Algorith)
    TIMESTRONG(ii)=timestrong; % Time taken for verifyimg strong duality property
    fprintf('test # %i and: Time Timestrong KKT1 KKT2 values are:  %g %g  %g %g \n',...
        ii,TIME(ii),TIMESTRONG(ii),KKT1(ii),KKT2(ii))
    [fval2,time2]=socpsdo(A,a,b,beta,delta,n);
    fbest=min([fval fval2]);
    ACCURACYNEW(ii)=(fval-fbest)/abs(fbest);  % accuracy of new algoeithm (Main Algorithm)
    ACCURACYSOCP(ii)=(fval2-fbest)/abs(fbest);  % accuracy of SOCP/SDO reformulation
    TIMESOCP(ii)=time2;    % time taken for solving SOCP/SDO reformulation
   
end
%====================================================================
% =========    averaged results  ====================================

TIMEs(nn)=sum(TIME)/notests;
TIMESTRONGs(nn)=sum(TIMESTRONG)/notests;
KKT1s(nn)=sum(KKT1)/notests;
KKT2s(nn)=sum(KKT1)/notests;
ACCURACYNEWs(nn)=sum(ACCURACYNEW)/notests;
ACCURACYSOCPs(nn)=sum(ACCURACYSOCP)/notests;
TIMESOCPs(nn)=sum(TIMESOCP)/notests;
nLNGMs(nn)=sum(nLNGM);    % number of test problems for which the optimal solution is LNGM
fprintf('\n \n #tests/dim %i and average Time, Timestrong, KKT1, KKT2, accuracynew, accuracysocp values are:  %g %g  %g %g %g %g\n',...
    notests,mean(TIME),mean(TIMESTRONG),mean(KKT1),mean(KKT2),mean(ACCURACYNEW),mean(ACCURACYSOCP))
fprintf('number of tests %i and number of NLGM optimal %i \n',notests,nLNGM)

end   % for nn = loop
outputdatamatrixT2=[ TIMEs' TIMESTRONGs' KKT1s' KKT2s' ACCURACYNEWs' ACCURACYSOCPs' TIMESOCPs' nLNGMs' ];
save('matrixoutputT2','outputdatamatrixT2','ns')
TableResultsT2;

% profile report
