%%%% TestScriptPrecise
%
%Script to generate the plots presented in Section 5 of the manuscript
%"Error Bounds and Singularity Degree in Semidefinite Programming."
%Each memeber of the matrix sequences X(sigma^k) and Z(sigma^k) is
%obtained by solving an SDP.  While this approach is more accurate, it is
%much more time consuming than the approach of TestScript.  Moreover,
%we have not observed an advantage to this approach as far detecting poor
%convergence and identifying singularity degree.
%
%Sremac, Woerdeman, and Wolkowicz
%2020

clear
seed=1;
rng(seed);
%rng('shuffle');         % comment out if same problems desired
rngsave = rng;          % comment out if same problems desired
save('rngfile','rngsave')
fprintf('\n Generate the SECOND set of plots in Section 5 of the manuscript\n')
fprintf('For each memeber of matrix sequences X(sigma^k) and Z(sigma^k)\n')


n = 10;             %set dimension parameters
tn = n*(n+1)/2;     %triangular number
tol = 12;           %stopping criteria tolerance (in decimal spots)
ptype = 'WNM';     %TUN,TOE,WNM,WEWO


%% for ssMat/ssvec
n2=n^2;
svind = find(triu(ones(n),1)); %linear indices of the strictly upper triangle
sMind = zeros(n2,1); %matrix indices with respect to linear indices of the diagonal and the strictly upper triangle
sMind(svind) = 1:tn-n; %strictly upper triangle
sMind(1:n+1:n2) = tn-n+1:tn;%the diagonal

%calculate the linear indices of the strictly lower triangle
indM = zeros(n); %matrix indices
indM(:) = 1:n2;
indM = indM';
lsvind = indM(svind); %strictly lower triangular
sMind(lsvind) = 1:tn-n; %strictly lower triangular

%% Problem generation
if strcmp(ptype,'TUN')==1     %%%%%%% Generalized Tuncel Problem
    
    Rot = eye(n);       %Orthogonal matrix, in case isomorphism desired
    SDTun = 4;       %singularity degree desired
    singdeg = SDTun;
    tracebound = 1;     %bound for trace to get bdd problem
    [A,b,Xorig] = MakeProblemTun(n,tracebound,svind,SDTun,Rot);
    [m,~] = size(A);
    AA = cell(m,1);
    for i=1:m
        AA{i} = ssMat(A(i,:)',n,sMind);
    end
    
elseif strcmp(ptype,'TOE')==1     %%%%%%%% Toeplitz Cycle Completion
    
    [A,b,Xorig,Exp,theta,phi] = MakeProblemToe(n,2,svind,pi/9,pi);
    m = size(A,1);
    
elseif strcmp(ptype,'WNM')==1    %%%%%%%% Waki-Muramatsu Problem
    
    singdeg = 5;
    r = (n-1)/3;
    if r>floor(r)
        fprintf('Poor choice of n for problem WNM \n');
        keyboard
    end
    [A,b,Xorig] = MakeProblemWNM(r,svind);
    m = size(A,1);
    AA = cell(m,1);
    for ii=1:m
        AA{ii} = ssMat(A(ii,:),n,sMind);
    end
    
elseif strcmp(ptype,'WEWO')==1    %%%%%%%% Wei-Wolkowicz Problem Generation
    
    dualHasInteriorFlag=1;      %for dual strict feasibility
    m = 2*n;
    NoP=round(n/3);             %maximal primal rank
    gap=6;                      %complementarity gap
    NoD=n-NoP-gap;              %maximal dual rank
    [AA, bb, C, Xorig, yorig, Zorig, status] = ...
        genBlockHard(m-1, n, gap, NoP, NoD, dualHasInteriorFlag);
    A = zeros(m,tn);
    for i=1:m-1
        A(i,:) = ssvec(AA{i},svind);
        AA{i} = ssMat(A(i,:)',n,sMind);
    end
    A(m,:) = ssvec(C,svind)';
    AA{m} = C;
    b = A*ssvec(Xorig,svind);
    
end

%% loop to create matrix sequences X(sigma^k) and Z(sigma^k)
sigma = 0.6;
sigmavec = zeros(60,1);
for ii=1:length(sigmavec)   % linear decrease
    sigmavec(ii) = sigma^(ii-1);
end

eigXstore = zeros(length(sigmavec),n);
eigZstore = eigXstore;
diagZstore = eigZstore;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(sigmavec)   % main for loop for decreasing alpha
    
    data.A=A;
    data.b=b + sigmavec(i)*A*ssvec(eye(n),svind);
    data.m=m;
    data.n=n;
    data.X0=Xorig + sigmavec(i)*eye(n);
    params.quiet=true;
    params.maxiter=63;
    params.maxlastiter=-1;              %set to negative if no extra iters desired
    params.svind=svind;
    params.sMind=sMind;
    params.tolalpha = 10^(-tol);        %tolerance for alpha
    params.tolResNorm = 10^(-tol);      %tolerance for residuals
    
    [output,psdflag] = IPGN(data,params);
    eigXstore(i,:) = eig(output.X)';
    eigZstore(i,:) = eig(output.Z)';
    diagZstore(i,:) = diag(output.Z)';
end

%% Plots of Eigenvalue Ratios
%R_Q(i,k) approximate plot
forplot1 = eigXstore(2:end,:)./eigXstore(1:end-1,:);
figure(1)
clf
hold on
plot(ones(length(eigXstore)-1,1),...
    'Color',[0.5 0.5 0.5],'LineStyle','-')
for kk=1:n-1
    plot(sigma^(2^(-(kk-1)))*ones(length(eigXstore)-1,1),...
        'Color',[0.5 0.5 0.5],'LineStyle','--')
end
for jj=1:n
    plot(forplot1(:,jj),'-')
end
ax = gca;
ax.FontSize = 16;
title('\textbf{Plot of $R_Q(i,k)$}',...
    'interpreter','latex','FontSize',20)
xlabel('$k$','interpreter','latex')
ylabel('$R_Q(i,k)$',...
    'interpreter','latex')
hold off

%R_N(i,k) approximate plot
forplot2=zeros(size(eigXstore,1)-1,size(eigXstore,2)-1);
for ii=2:size(eigXstore,2)
    forplot2(:,ii-1)=eigXstore(2:end,ii)./eigXstore(2:end,ii-1);
end
figure(2)
clf
hold on
for ii=1:size(forplot2,2)
    plot(log10(forplot2(:,ii)),'-')
end
ax = gca;
ax.FontSize = 16;
title('\textbf{Plot of $R_N(i,k)$}','interpreter','latex',...
    'FontSize',20)
xlabel('$k$','interpreter','latex')
ylabel('$\log_{10}(R_N(i,k))$',...
    'interpreter','latex')
hold off

%Convergence rates of diag(Z)
forplot3=zeros(size(diagZstore,1)-1,size(diagZstore,2)-1);
for ii=2:size(diagZstore,2)
    forplot3(:,ii-1)=diagZstore(2:end,ii)./diagZstore(2:end,ii-1);
end
figure(3)
clf
hold on
for ii=1:size(forplot3,2)
    plot(log10(forplot3(:,ii)),'-')
end
ax = gca;
ax.FontSize = 16;
title('\textbf{Plot of diag$(Z(\sigma^k))$}','interpreter','latex',...
    'FontSize',20)
xlabel('$k$','interpreter','latex')
ylabel('diag$(Z(\sigma^k))$',...
    'interpreter','latex')
hold off

%Q-convergence rates of diag(Z)
forplot4 = diagZstore(2:end,:)./diagZstore(1:end-1,:);
figure(4)
clf
hold on
plot(ones(length(diagZstore)-1,1),...
    'Color',[0.5 0.5 0.5],'LineStyle','-')
for kk=1:n-1
    plot(sigma^(2^(-(kk-1)))*ones(length(eigXstore)-1,1),...
        'Color',[0.5 0.5 0.5],'LineStyle','--')
end
for jj=1:n
    plot(forplot1(:,jj),'-')
end
ax = gca;
ax.FontSize = 16;
title('\textbf{Plot of $Q$-convergence of diag$(Z(\sigma^k))$}',...
    'interpreter','latex','FontSize',20)
xlabel('$k$','interpreter','latex')
ylabel('diag$(Z(\sigma^k))$',...
    'interpreter','latex')
hold off
