% This codes use primal interior point method to solve the optimal design 
% with D criterion, namely, 
% 
% minimize {logdet(M^(-1)): sum(w) = 1, w >= 0, M(w)>0}, 
% where M(w) = w_1*A_1 + ... + w_m*A_m, A_i's are n by n psd matrices.
%%
% Input:
% M        -  a matrix whose ith column consists of the vector formed by 
%             stacking the columns of A_i one by one; 
% tol      -  tolerance parameter;
% maxiter  -  maximum number of iterations allowed.
%
%%
% Output:
% w     - approximate optimal solution
% fval  - approximate optimal value
% iter  - number of iterations
% time  - cpu time
%
%%
% 
% Written by Zhaosong Lu and Ting Kei Pong
% Jan 12, 2012

function [w,fval,iter,time] = IPM_D(M,tol,maxiter)

tic

% initialization
mu_final = tol;
m = size(M,2);
n = sqrt(size(M,1));
w = ones(m-1,1)/m;
M = rearrange(M);
[MA,c,index] = selectw(M);
blk{1,1} = 's';
blk{1,2} = n;

ME = [MA; ones(1,m-1)];

mu = 10;
iter = 1;
dw = zeros(m-1,1);

while mu >= mu_final && iter < maxiter
  
  X = smat(blk,MA*w + c);
  X = (X + X')/2;
  
  fval = function_value(X,w,mu);  
  [g,invX] = gradient(X,w,MA,mu,blk);  
  g_tol = max(tol,mu);
  normg = norm(g);
  
  while normg > g_tol && iter < maxiter
    
    dw = search_direction(X,blk,w,MA,ME,mu,g,invX,dw);
    dX = smat(blk,MA*dw);
    
    stepsize = 1;
    mindw = min(dw./w);
    if mindw < 0
      stepsize = min([stepsize,-0.95/mindw,1]);
    end
    sumdw = sum(dw);
    if sumdw > 0
      stepsize = min([stepsize,0.95*(1-sum(w))/(sumdw),1]);
    end

    Xnew = X + stepsize*dX;
    wnew = w + stepsize*dw;
    Xnew = (Xnew + Xnew')/2;
    
    fnew = function_value(Xnew,wnew,mu);    
    dderiv = g'*dw;
    
    while fnew >= fval + 0.1*stepsize*dderiv && stepsize > 1e-8 && dderiv < 0 % Armijo rule
      stepsize = stepsize/2;
      
      Xnew = X + stepsize*dX;
      wnew = w + stepsize*dw;
      Xnew = (Xnew + Xnew')/2;
      
      fnew = function_value(Xnew,wnew,mu);
    end
    
    if stepsize <= 1e-8 || dderiv >= 0
%       fprintf(' small stepsize size\n')
      break
    end
    
    % Update the variables
    
    X = Xnew;
    w = wnew;
    
    if abs(fval - fnew)/max(abs(fval),1) < min(tol,1e-10)
      break
    end
    
    fval = fnew;
    [g,invX] = gradient(X,w,MA,mu,blk);
    normg = norm(g);
    
    % fprintf(' mu, iter, fval, dderiv, stepsize %g, %g, %g, %g, %g\n', mu, iter, fval, dderiv, stepsize)
    
    iter = iter + 1;
    
  end
  
%   fprintf(' mu: %5.1e    iter: %d    fval: %g    stepsize: %5.1e\n', mu, iter, fval, stepsize)
  
  mu = mu/2;
  
end

wnew = zeros(m,1);
wnew(1:index-1) = w(1:index-1);
wnew(index) = 1 - sum(w);
wnew(index+1:m) = w(index:m-1);

w = wnew;

D = svd(smat(blk,M*w));
fval = -sum(log(D));

fprintf(' iter: %d   optimal_value: %g   \n', iter, fval)

time = toc;
    