% This code computes a matrix in the neighborhood (not containing zero)
% (Q_ij + eta_ij)>= X_ij >= (Q_ij - eta_ij)
% with minimum condition number, where Q is a positive definite
% n by n matrix and eta is symmetrix. It applies the alternating
% direction method to solve
%
% min  lambda_max(X)
% s.t. X - I >= 0, t >= 0, (Q_ij + eta_ij)t>= X_ij >= (Q_ij - eta_ij)t
%% 
% Input
% Q       - a positive definite matrix
% eta     - a symmetric matrix so that the neighborhood does not contain
%           the zero matrix
% tol     - tolerance parameter
% maxiter - maximum number of iterations allowed
%%
% Output
% answer  - approximate minimum condition number
% X       - matrix in the convex hull with cond(X) = t
% time    - cpu time
% iter    - number of iterations
%%
% Subroutines called:
% update_X_adm_bd, update_Yt_adm_bd, fval_cond_bd, grad_cond_bd, simplex_proj
%%
% March 5, 2011
% Written by Zhaosong Lu and Ting Kei Pong

function [answer,X,iter,time] = cond_adm_bd(Q,eta,tol,maxiter)

tic
% parameters
n = size(Q,1);
beta = 0.1;

% initialization
Y = eye(n);
P = -eye(n)/n;
w = ones(n,1)/n/beta;
t = 1;
lb = Q - eta;
ub = Q + eta;

iter = 1;

while iter <= maxiter
  
  [X,w] = update_X_adm_bd(n,w,Y,P,beta,tol/100);
  X = (X + X')/2;
  
  [Y,t] = update_Yt_adm_bd(n,lb,ub,X,Y,P,beta,t,tol/100);
  Y = (Y + Y')/2;
  
  P = P + beta*(X - Y);
  P = (P + P')/2;
  
  eigY = eig(Y);
  pval = max(eigY);
  pfeas = (1 - min(eigY))/max(1,norm(eigY)); % violations of interval constraints
  
  p = eig(P);
  dval = 1 + sum(p);
  dfeas1 = (sum(max(-p,0)) - 1)/max(1,norm(p));
  dfeas2 = (sum(sum((min(P,0)).*lb)) + sum(sum((max(P,0)).*ub)))/max(1,norm(p));
  
  if abs(pval - dval)/max(1,abs(pval)) < 10*tol && max([dfeas1,dfeas2]) < tol && pfeas < tol && min(eigY) > 0
    fprintf('TERMINATE:\n iter = %6d, pval = %6.1f, dval = %6.1f, pfeas = %6.5f, dfeas1 = %6.5f, dfeas2 = %6.5f, time = %4.1f\n',...
      iter, pval, dval, pfeas, dfeas1, dfeas2, toc)
    break
  end
  
%   if mod(iter,2000) == 0
%     fprintf(' iter = %6d, pval = %6.1f, dval = %6.1f, pfeas = %6.5f, dfeas1 = %6.5f, dfeas2 = %6.5f, time = %4.1f\n',...
%       iter, pval, dval, pfeas, dfeas1, dfeas2, toc)
%   end
  
  iter = iter + 1;
 
end

X = Y/t;
answer = cond(X);
time = toc;
fprintf(' condition number = %4.1f, time = %4.1f\n', answer, time)