% This code constructs input to SDPT3 for A-criterion:
%
% minimize {tr(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 matrices and M^+ is the pseudoinverse of M.
% 
% Input
% M - matrix whose ith column is formed by 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] = A_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} = n + k;
block{1} = 's'; block{2} = n;

nk = (n+k)*(n+k+1)/2;

Aw = zeros(n*(n+1)/2 + n*k,m);
for i = 1:m
  tmp = reshape(M(:,i),n,n);
  Mw = svec(block,tmp);
  Aw(:,i) = [zeros(n*k,1); Mw];
end
Aw = [Aw;ones(1,m)];

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

AX1 = zeros(n*k,nk);
AX2 = zeros(n*(n+1)/2,nk);

count = 1;
for j = 1:n
  for i = 1:k
    E = zeros(k,n);
    E(i,j) = 1/2;
    tmp = [zeros(k,k) E;zeros(n,k) zeros(n,n)];
    AX1(count,:) = svec(blk(2,:),tmp)';
    count = count + 1;
  end
end

count = 1;
for j = 1:n
  for i = 1:j
    E = zeros(n,n);
    if i ~= j
      E(i,j) = 1/sqrt(2);
      E(j,i) = 1/sqrt(2);
    else E(j,j) = 1;
    end
    tmp = [zeros(k,k) zeros(k,n);zeros(n,k) E];
    AX2(count,:) = svec(blk(2,:),tmp)';
    count = count + 1;
  end
end

AX = [AX1; -AX2; zeros(1,nk)];

A{2,1} = AX';
C{2,1} = [eye(k) zeros(k,n); zeros(n,n+k)];

b = [reshape(K',n*k,1);svec(block,zeros(n,n));1];

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

% Y = X;

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;
D = svd(KtinvX*K);

fval = sum(D);

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