% runcode
clc
path('./gp_codes',path)
clear
n = 250;

Lambda = n*[10 100 1000];
robust = [0.05 0.1 0.5 1]; % the closer to being just infeasible, the slower.
fid = fopen('real_datap_250.txt','w');

% load given data matrix
load cov6
G0 = cov6;
G0 = (G0 + G0')/2;
% G0 = G0(1:n,1:n);
diagD = sqrt(diag(G0));
invdiagD = 1./diagD;
G0 = diag(invdiagD)*G0*diag(invdiagD); % scale to a correlation matrix
eigG0 = eig(G0);

n = size(G0,1);
k = ceil(n/2); % number of stocks movable
r = 1; % number of inequality constraints
m = ceil(n*(n-1)/40); % number of equality constraints, 2% of all entries
thresold = 0;%ceil(m/2);

% adding a weight representing the VaR
w = ones(n,1)/n;
Dw = diagD.*w;
D = diag(diagD);
DwwD = Dw*Dw';

fprintf('\n done generation part 1: min eig = %g, max eig = %g,\n original VaR = %g\n', min(eigG0), max(eigG0), norminv(0.995)*sqrt(sum(sum(G0.*DwwD))))

% generate random indices for equality constraints
Bindex = zeros(m,2);
Index = randperm(n*(n-1)/2);
Index = Index(1:m);
d = zeros(m,1);
for ii = 1:m
  tmp = Index(ii);
  [i,j] = sub2indsym(tmp);
  Bindex(ii,:) = [n*(i-1)+j n*(j-1)+i];
  if i == j
    fprintf(' error\n')
    return
  end
  if ii > thresold
    d(ii) = 0.9;
  else
    d(ii) = G0(i,j);
  end
end

% Generate random uncertainty centers
Index = randperm(n);
Index = Index(1:k); % number of centers = k - 1; number of stocks moved = k
Q = null(ones(1,k)); % basis for sum(w) = 0, k by k-1 matrix
M = zeros(n,k-1);
M(Index,:) = Q;
Mx = kron(Dw,D*M); % the multiple by 2 is taken into account in the symmetrization below
P = cell(1,1);
q = zeros(1,k);

fprintf(' done generation\n')
sol = cell(length(robust),length(Lambda));

for ii = 1:length(Lambda)
  q(1,1) = -sum(sum(DwwD.*G0));
  absq = abs(q(1,1)); % normalize data
  q(1,1) = sign(q(1,1));
  lambda = Lambda(ii); % heuristics
  G = G0 + lambda*DwwD;
  
  for jj = 1:length(robust)
    
    P{1} = -DwwD(:)'/absq;
    ek = robust(jj); % robust radius
    for j = 1:k-1
      A = -reshape(Mx(:,j),n,n)/absq;
      A = A + A';
      P{1} = [P{1};ek*A(:)'];
    end
    
    [fval, dbest, pfeas, gap, X, iter, time, exitflag] = lssdp_gp(G,P,q,Bindex,d,1e-9,25000); 
    pfeas = max(pfeas,-min(eig(X)));

    fprintf(fid,' %8d & %4.3f &', lambda, robust(jj));
    if exitflag == 1
      fprintf(fid,' %5.1f/%8d/%5.4e/%2.1e/%6.5e/%6.5e &', time, iter, fval, pfeas, norm(X - G0,'fro'), norminv(0.995)*sqrt(sum(sum(X.*DwwD))));
    elseif exitflag == -1
      fprintf(fid,' Infeasible: cpu/iter = %5.1f/%8d&', time, iter);
    elseif exitflag == 0
      fprintf(fid,' Early term: cpu/iter/fval/pfeas = %5.1f/%8d/%5.4e/%2.1e&', time, iter, fval, pfeas);
    end
      
      
    [fval, X, time] = lssdp_sdpt3_p(G,P,q,Bindex,d);
    
    pfeas = 0;
    % routine for computing gradient, taken from gradient.m
    B2 = Bindex(:,1);
    vecX = X(:);
    grad(1:k,1) = P{1}*vecX - q(1,:)';
    grad(k+1:k+n,1) = diag(X) - 1;
    grad(k+n+1:k+n+m,1) = 2*vecX(B2) - 2*d;
    grad = -grad;
      
    v = grad(1:k,1);
    pfeas = max(pfeas, norm(v - soc_proj(v)));
    pfeas = max([pfeas, norm(grad(k+1:k+n)), norm(grad(k+n+1:end))/2]); % absolute primal feasibility
    pfeas = max(pfeas,-min(eig(X)));
      
    if fval >= 0 
      fprintf(fid,' %5.1f/%5.4e/%2.1e/%6.5e/%6.5e\n', time, fval, pfeas, norm(X - G0,'fro'), norminv(0.995)*sqrt(sum(sum(X.*DwwD))));
    elseif fval < 0
      fprintf(fid,' Infeasible: cpu = %5.1f', time);
      time = lssdp_sdpt3_p_feas(G,P,q,Bindex,d);
      fprintf(fid,' (%5.1f)\n', time);
    end
  
  end
  
  
end
fclose(fid);


