function [x,E,z] = ipm4mve(A,b,x0)

%  Find the maximum volume ellipsoid
%     Ell = {v:  v = x + Es, ||s|| <= 1}
%  or Ell = {v:  ||E^{-1}(v-x)|| <= 1}
%  inscribing a full-dimensional polytope
%          {v:  Av <= b}
%  Input:  A, b --- defining the polytope
%   (Optinal x0 --- interior point, A*x0 < b)
%  Output:  x --- center of the ellipsoid
%           E --- matrix defining ellipsoid: {x + Es | ||s|| <= 1}
%                 (the output is a triangular matrix)
%           z --- constraint slacks

%--------------------------------------
% Yin Zhang, Rice University, 07/29/02
% Last modified: 10/15/16
% Last modified: 09/15/19
% Name  changed: 03/05/20
% Last modified: 10/09/20 (added "loud") 
%--------------------------------------

global loud

maxiter = 80; tol1 = 1.e-4; tol2 = 1.e-7;
[m, n] = size(A); 
%t0 = tic;

if isempty(loud), loud = false; end

if nargin < 3
  mve_chkdata(A,b);
  tic, [msg,x0,t] = mve_presolve(A,b,maxiter,tol1);
  if loud, fprintf('  End of Presolve ... t = %f\n',t); toc, end
  if msg(1) ~= 's', disp(msg); x = []; E = []; return; end
  %A=bsxfun(@times,1./(b-A*x0),A); b=ones(m,1); x0=zeros(n,1);
  [x,E,~,~,z] = mve_solver(A,b,x0,maxiter,tol2);
else
  mve_chkdata(A,b,x0);
  [x,E,~,~,z] = mve_solver(A,b,x0,maxiter,tol2);
end
%fprintf('\n\t[m, n] = [%i, %i]\n', m, n);
%fprintf('  Elapsed time: %g seconds\n', toc(t0));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mve_chkdata(A,b,x0)
% check consistency of data

%--------------------------------------
% Yin Zhang, Rice University, 07/29/02
% Last modified: 10/15/16
%--------------------------------------

 [rA,cA] = size(A); [rb,cb] = size(b);
 if cA < 2 
    error('A must have at least 2 columns'); 
 end
 if rA < cA+1
    error('A has too few rows'); 
 end
 if cb ~= 1 || rb ~= rA
    error('size of b mis-matches');
 end

 if nargin < 3, return; end

 [rx0,cx0] = size(x0);
 if cA ~= rx0 || cx0 ~= 1
    error('size of x0 mis-matches');
 end
 if any(b-A*x0) <= eps
    error('x0 is not interior');
 end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [msg,x,t,s,y] = mve_presolve(A,b,maxiter,tol)
% Solve LP: max t, s.t. Ax + t*e <= b

%--------------------------------------
% Yin Zhang, Rice University, 07/29/02
% Last modified: 10/15/16
%--------------------------------------
global loud

[m, n] = size(A);  bnrm = norm(b); 
o_m = zeros(m,1);  o_n = zeros(n,1);
e_m =  ones(m,1);  %e_n =  ones(n,1);

% initialize
x = o_n; y = e_m/m;
t = min(b) - 1; s = b - t;

dx = o_n; dxc = dx; %#ok<*NASGU>
ds = o_m; dsc = ds;
dy = o_m; dyc = dy;
dt = 0;   dtc = 0;
tau0 = 0.995; 
sigma0 = 0.2;
msg = 'successful';

p = 1:n;
if issparse(A)
    absA = abs(A); 
    p = symamd(absA'*absA);
end
p(n+1) = n+1;

if loud
    fprintf('\n  Residuals:   Primal     Dual     Duality    Obj\n');
    fprintf('  --------------------------------------------------\n');
end

for iter = 0:maxiter

  % KKT residuals
  r1 = b - (A*x + s + t);
  r2 = -A'*y;
  r3 = 1 - sum(y);
  r4 = -s.*y;
  gap = -sum(r4);

  % relative residual norms and gap
  prif = norm(r1)/(1 + bnrm); 
  drif = norm([r2;r3]); 
  rgap = abs(b'*y - t)/(1 + abs(t)); 
  total_err = max([prif drif rgap]);
  
  % progress output & check stopping
  if (total_err < tol)
      if loud
          fprintf('  iter. %3i: %9.1e %9.1e %9.1e %9.1e\n',...
              iter,prif,drif,rgap,t);
      end
      break;
  end
  if loud
      fprintf('  iter. %3i: %9.1e %9.1e %9.1e %9.1e\n',...
          iter,prif,drif,rgap,t);
  end
  if dt > 1.e+3*bnrm || t > 1.e+6*bnrm
      msg = 'unbounded?'; break; 
  end  

  % Shur complement matrix
  d = min(5.e+15,y./s); 
  AtD = A'*spdiags(d,0,m,m);
  AtDe = AtD*e_m;
  AtDA = AtD*A;
  B = [AtDA AtDe; AtDe' sum(d)];
  B = B + 1.e-14*speye(n+1);

  % Cholesky decomposition
  if ~issparse(A)
      R = chol(B);
  else
      %R = cholinc(B(p,p),'inf');
      R = chol(B(p,p));
  end
  
  % predictor step & length
  [dx,ds,dt,dy] = calcstep(A,R,p,s,y,r1,r2,r3,r4);   

  alphap = -1/min([-1; ds./s]);
  alphad = -1/min([-1; dy./y]);
  
  % determine mu
  ratio = (s+alphap*ds)'*(y+alphad*dy)/gap;
  sigma = min(sigma0, ratio^2);
  mu = sigma*gap/m;
  
  % corrector and combined step & length
  [dxc,dsc,dtc,dyc] = calcstep(A,R,p,s,y,o_m,o_n,0,mu-ds.*dy);  
  dx = dx + dxc; ds = ds + dsc; 
  dt = dt + dtc; dy = dy + dyc;
  alphap = -1/min([-.5; ds./s]);
  alphad = -1/min([-.5; dy./y]);
  
  % update iterates
  tau = max(tau0, 1-gap/m);
  alphap = min(1,tau*alphap);
  alphad = min(1,tau*alphad);
  x = x + alphap * dx;
  s = s + alphap * ds;
  t = t + alphap * dt; 
  y = y + alphad * dy;
end
if t < sqrt(eps), msg = 'no volume'; end
 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dx,ds,dt,dy] = calcstep(A,R,p,s,y,r1,r2,r3,r4)
 dxdt = zeros(size(R,1),1);
 tmp = (r1.*y-r4)./s;
 rhs = [r2+A'*tmp; r3+sum(tmp)]; 
 if ~issparse(A)
     dxdt = R\(R'\rhs);
 else
     dxdt(p) = R\(R'\rhs(p));
 end
 dx = dxdt(1:end-1); 
 dt = dxdt(end);
 ds = r1 - A*dx - dt;
 dy = (r4 - y.*ds)./s;
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,E,msg,y,z,iter] = mve_solver(A,b,x0,maxiter,tol)
%  Find the maximum volume ellipsoid
%    {v:  v = x + Es, ||s|| <= 1}
%  inscribing a full-dimensional polytope
%          {v:  Av <= b}
%  Input:  A, b --- defining the polytope
%          x0 --- interior point (Ax0 < b)
%  Output: x --- center of the ellipsoid
%          E2 --- E'*E

%--------------------------------------
% Yin Zhang, Rice University, 07/29/02
% Last modified: 10/15/16
%--------------------------------------
global loud

[m, n] = size(A);
Anrm = norm(A,'inf'); 
bnrm = norm(b,1); 

if ~exist('maxiter','var'), maxiter = 50; end
if ~exist('tol','var'), tol = 1.e-4; end
minmu = 2*eps; 
tau0 = .99; mu0 = 0.1;

bmAx0 = b - A*x0;
if any(bmAx0<=0), error('x0 not interior'); end

A = sparse(1:m,1:m,1./bmAx0)*A; b = ones(m,1); 
x = zeros(n,1); y = ones(m,1); bmAx = b;

if loud == true
    fprintf('\n  Residuals:   Dual     Primal    Duality  logdet(E)\n');
    fprintf('  --------------------------------------------------\n');
end

%res = 1; 
msg = 0;
for iter=1:maxiter %----- loop starts -----

if iter > 1, bmAx = bmAx - astep*Adx; end

YA = bsxfun(@times,y,A);
AtYA = A'*YA;
AtYA = full(AtYA);
% E2 = inv(AtYA); objval = log(det(E2))/2;
U = chol(AtYA); d = diag(U);
%objval = -log(prod(full(d)));
objval = -sum(log(d));
AUi = A / U;
% Q = A*E2*A'; 
Q = AUi * AUi';
h = sqrt(diag(Q));

if iter == 1
   t = min(bmAx./h); 
   y = y/t^2; h = t*h;
   z = max(1.e-1, bmAx-h);
   Q = t^2*Q; YA = YA/t^2;
end

%% residuals & stop rule
yz = y.*z; yh = y.*h;
gap = sum(yz)/m;
rmu = min(mu0, gap)*gap;
rmu = max(rmu, minmu);

R1 = -(yh'*A)';
R2 = bmAx - h - z;
R3 = rmu - yz;

r1 = norm(R1,'inf')/Anrm;
r2 = norm(R2,'inf')/bnrm;
r3 = gap;
%r3 = norm(R3,'inf');
res = max([r1 r2 r3]);

if loud
    fprintf('  iter %3i  ', iter);
    fprintf('%9.1e %9.1e %9.1e  %12.8e\n',r1,r2,r3,objval);
end
if res < tol && rmu <= minmu 
   if loud, fprintf('  Converged!\n'); end
   x = x + x0; msg = 1; break; 
end

%% compute steps
YQ = bsxfun(@times,y,Q); 
YQQY = YQ.*YQ'; yh2 = 2*yh;
M  = YQQY + sparse(1:m,1:m,max(yh2.*z,eps/8));
M = chol(M); 
YA = bsxfun(@times,y,A);
T = M \ (M' \ bsxfun(@times,h+z,YA));
AtNMi = (bsxfun(@times,yh2,T)-YA)';
AtNMiA = AtNMi*A;% + sqrt(2*eps)*eye(n);
R3Dy = R3./y; R23 = R2 - R3Dy;
warning('off', 'MATLAB:nearlySingularMatrix')
dx = AtNMiA \ (R1 + AtNMi*R23); 
Adx = A*dx;
yidy = M \ (M' \ (yh2.*(Adx - R23)));
dy = y.*yidy;
dz = R3Dy - z.*yidy;

%% step length & update
ax = -1/min([-Adx./bmAx; -.5]);
ay = -1/min([ yidy; -.5]);
az = -1/min([dz./z; -.5]); 
tau = max(tau0, 1-gap);
astep = min(1,tau*min([ax ay az]));

x = x + astep*dx;
y = y + astep*dy;
z = z + astep*dz;

end
% E2 = inv(AtYA); U = chol(AtYA);
E = inv(U);
end