% This code constructs input to SDPT3 for D-criterion:
%
% minimize { logdet(K'*(M(w))^(+)*K) : e^Tw=1, w>=0, Range(K) \subseteq Range(M)},
% where K is an n by k matrix, M(w) = w_1A_1+...+w_mA_m,
% A_i are n by n psd matrices, and M^+ is the pseudoinverse of M.
% 
% Input
% M - matrix whose ith column is formed from stacking columns of A_i together;
% K - n by k matrix;
%
% Output
% w    - approximate optimal solution;
% fval - approximate optimal value;
% time - time for running SDPT3
%
% Jan 12, 2012
% Written by Zhaosong Lu and Ting Kei Pong

function [w, fval, time] = D_criterion_SDPT3_K(M,K)

n = sqrt(size(M,1));
m = size(M,2);
k = size(K,2);

blk{1,1} = 'l'; blk{1,2} = m;
blk{2,1} = 's'; blk{2,2} = k;
blk{3,1} = 's'; blk{3,2} = n;

count = 1;
MA = zeros(n*(n+1)/2,m);
for j = 1:n
  for i = 1:j
    if i == j
      MA(count,:) = M(n*(i-1) + j,:);
      count = count + 1;
    end
    if i < j
      MA(count,:) = sqrt(2)*M(n*(i-1) + j,:);
      count = count + 1;
    end
  end
end

Aw = [MA;ones(1,m)];

A{1,1} = Aw';
C{1,1} = zeros(m,1);

Au = zeros(n*(n+1)/2,k*(k+1)/2);
count = 1;
for j = 1:k
  for i = 1:j
    E = zeros(k,k);
    if i ~= j
      E(i,j) = 1/sqrt(2);
      E(j,i) = 1/sqrt(2);
    else E(j,j) = 1;
    end
    KUKT = K*E*K';
    Au(:,count) = svec(blk(3,:),-KUKT);
    count = count + 1;
  end
end
Au = [Au;zeros(1,k*(k+1)/2)];

A{2,1} = Au';
C{2,1} = zeros(k,k);

A{3,1} = [-eye(n*(n+1)/2);zeros(1,n*(n+1)/2)]';
C{3,1} = zeros(n,n);

b = [zeros(n*(n+1)/2,1);1];

OPTIONS.parbarrier{1} = [];
OPTIONS.parbarrier{2} = 1;
OPTIONS.parbarrier{3} = [];

OPTIONS.printlevel = 0;
tic
[obj,X] = sdpt3(blk,A,C,b,OPTIONS); 
time = toc;

w = X{1};
w = simplex_proj(w,1);

X = reshape(M*w,n,n);

[V,D] = svd(X);
invX = V*diag(diag(D).^(-1))*V';
KtinvX = K'*invX;
W = KtinvX*K;
D = svd(W);

fval = sum(log(D));   

fprintf(' fval = %g   \n', fval)