% Boyd & Vandenberghe, "Convex Optimization"
% Original version by Lieven Vandenberghe
% Updated for CVX by Almir Mutapcic - Jan 2006
% (a figure is generated)
%
% This is an example of D-optimal, A-optimal, and E-optimal
% experiment designs.

% problem data
m = 10;
angles1 = linspace(3*pi/4,pi,m);
angles2 = linspace(0,-pi/2,m);

% sensor positions
V = [3.0*[cos(angles1); sin(angles1)], ...
     1.5*[cos(angles2); sin(angles2)]];
p = size(V,2);
n = 2;
noangles = 5000;

% D-optimal design
%
%      maximize    log det V*diag(lambda)*V'
%      subject to  sum(lambda)=1,  lambda >=0
%

% setup the problem and solve it
cvx_begin
  variable lambda(p)
  maximize ( det_rootn( V*diag(lambda)*V' ) )
  subject to
    sum(lambda) == 1;
    lambda >= 0;
cvx_end
lambdaD = lambda; % save the solution for confidence ellipsoids

% plot results
figure(1)
% draw ellipsoid v'*W*v <= 2
W = inv(V*diag(lambda)*V');
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(2)*(R\[cos(angles); sin(angles)]);
d = plot(ellipsoid(1,:), ellipsoid(2,:), '--', 0,0,'+');
set(d, 'Color', [0 0.5 0]); set(d(2),'MarkerFaceColor',[0 0.5 0]);
hold on;

dot=plot(V(1,:),V(2,:),'o');
ind = find(lambda > 0.001);
dots = plot(V(1,ind),V(2,ind),'o');
set(dots,'MarkerFaceColor','blue');

% print out nonzero lambda
disp('Nonzero lambda values for D design:');
for i=1:length(ind)
   text(V(1,ind(i)),V(2,ind(i)), ['l',int2str(ind(i))]);
   disp(['lambda(',int2str(ind(i)),') = ', num2str(lambda(ind(i)))]);
end;

%axis([-4.5 4.5 -4.5 4.5])
axis([-5 5 -5 5])
set(gca,'Xtick',[]);
set(gca,'Ytick',[]);
hold off, axis off
% print -deps Ddesign.eps

% A-optimal design
%
%      minimize    Trace (sum_i lambdai*vi*vi')^{-1}
%      subject to  lambda >= 0, 1'*lambda = 1
%

% SDP formulation
e = eye(2,2);
cvx_begin sdp
  variables lambda(p) u(n)
  minimize ( sum(u) )
  subject to
    for k = 1:n
      [ V*diag(lambda)*V'  e(:,k);
        e(k,:)             u(k)   ] >= 0;
    end
    sum(lambda) == 1;
    lambda >= 0;
cvx_end
lambdaA = lambda; % save the solution for confidence ellipsoids

% plot results
figure(2)
% draw ellipsoid v'*W*v <= mu
W = inv(V*diag(lambda)*V')^2;
mu = diag(V'*W*V);
mu = mean(mu(ind));
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(mu)*(R\[cos(angles); sin(angles)]);
d = plot(ellipsoid(1,:), ellipsoid(2,:), '--',0,0,'+');
set(d, 'Color', [0 0.5 0]);
set(d(2), 'MarkerFaceColor', [0 0.5 0]);
hold on

dot = plot(V(1,:),V(2,:),'o');
ind = find(lambda > 0.001);
dots = plot(V(1,ind),V(2,ind),'o');
set(dots,'MarkerFaceColor','blue');

disp('Nonzero lambda values for A design:');
for i=1:length(ind)
   text(V(1,ind(i)),V(2,ind(i)), ['l',int2str(ind(i))]);
   disp(['lambda(',int2str(ind(i)),') = ', num2str(lambda(ind(i)))]);
end;
%axis([-4.5 4.5 -4.5 4.5])
axis([-5 5 -5 5])
set(gca,'Xtick',[]);
set(gca,'Ytick',[]);
axis off, hold off
% print -deps Adesign.eps

% E-optimal design
%
%      minimize    w
%      subject to  sum_i lambda_i*vi*vi' >= w*I
%                  lambda >= 0,  1'*lambda = 1;
%

cvx_begin sdp
  variables t lambda(p)
  maximize ( t )
  subject to
    V*diag(lambda)*V' >= t*eye(n,n);
    sum(lambda) == 1;
    lambda >= 0;
cvx_end

lambdaE = lambda; % save the solution for confidence ellipsoids

figure(3)
% draw ellipsoid v'*W*v <= mu
mu = diag(V'*W*V);
mu = mean(mu(ind));
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(mu)*(R\[cos(angles); sin(angles)]);
d = plot(ellipsoid(1,:), ellipsoid(2,:), '--', 0, 0, '+');
set(d, 'Color', [0 0.5 0]);
set(d(2), 'MarkerFaceColor', [0 0.5 0]);
hold on

dot = plot(V(1,:),V(2,:),'o');
lambda = lambda(1:p);
ind = find(lambda > 0.001);
dots = plot(V(1,ind),V(2,ind),'o');
set(dots,'MarkerFaceColor','blue');

disp('Nonzero lambda values for E design:');
for i=1:length(ind)
   text(V(1,ind(i)),V(2,ind(i)), ['l',int2str(ind(i))]);
   disp(['lambda(',int2str(ind(i)),') = ', num2str(lambda(ind(i)))]);
end;
%axis([-4.5 4.5 -4.5 4.5])
axis([-5 5 -5 5])
set(gca,'Xtick',[]);
set(gca,'Ytick',[]);
axis off, hold off
% print -deps Edesign.eps


% confidence ellipsoids
eta = 6.2514; % chi2inv(.9,3) value (command available in stat toolbox)
% draw 90 percent confidence ellipsoid  for D design
W = V*diag(lambdaD)*V';
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(eta)*(R\[cos(angles); sin(angles)]);

figure(4)
plot(0,0,'ok',ellipsoid(1,:), ellipsoid(2,:), '-');
text(ellipsoid(1,1100),ellipsoid(2,1100),'D');
hold on

% draw 90 percent confidence ellipsoid  for A design
W = V*diag(lambdaA)*V';
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(eta)*(R\[cos(angles); sin(angles)]);
plot(0,0,'ok',ellipsoid(1,:), ellipsoid(2,:), '-');
text(ellipsoid(1,1),ellipsoid(2,1),'A');

% draw 90 percent confidence ellipsoid  for E design
W = V*diag(lambdaE)*V';
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(eta)*(R\[cos(angles); sin(angles)]);
d=plot(0,0,'ok',ellipsoid(1,:), ellipsoid(2,:), '-');
set(d,'Color',[0 0.5 0]);
text(ellipsoid(1,4000),ellipsoid(2,4000),'E');

% draw 90 percent confidence ellipsoid  for uniform design
W_u = inv(V*V'/p);
R = chol(W_u);  % W = R'*R
ellipsoid_u = sqrt(eta)*(R\[cos(angles); sin(angles)]);
plot(ellipsoid_u(1,:), ellipsoid_u(2,:), '--');
text(ellipsoid_u(1),ellipsoid_u(2),'U');
set(gca,'Xtick',[]);
set(gca,'Ytick',[]);
axis off
% print -deps confidence.eps
hold off
 
Calling sedumi: 33 variables, 10 equality constraints
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 10, order n = 27, dim = 41, blocks = 3
nnz(A) = 91 + 0, nnz(ADA) = 88, nnz(L) = 49
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            2.75E+02 0.000
  1 :  -8.93E-01 7.75E+01 0.000 0.2814 0.9000 0.9000   1.64  1  1  3.2E+01
  2 :  -1.58E+00 2.57E+01 0.000 0.3321 0.9000 0.9000   1.17  1  1  1.4E+01
  3 :  -2.52E+00 9.03E+00 0.000 0.3509 0.9000 0.9000   0.61  1  1  5.8E+00
  4 :  -2.81E+00 3.10E+00 0.000 0.3436 0.9000 0.9000   1.01  1  1  2.0E+00
  5 :  -3.09E+00 7.85E-01 0.000 0.2531 0.9000 0.9000   0.84  1  1  5.4E-01
  6 :  -3.18E+00 3.50E-02 0.000 0.0446 0.9900 0.9900   0.98  1  1  2.4E-02
  7 :  -3.18E+00 8.11E-04 0.000 0.0232 0.9900 0.9900   1.00  1  1  5.6E-04
  8 :  -3.18E+00 2.94E-05 0.000 0.0362 0.9900 0.9900   1.00  1  1  2.0E-05
  9 :  -3.18E+00 2.99E-06 0.000 0.1019 0.9064 0.9000   1.00  1  1  2.4E-06
 10 :  -3.18E+00 3.37E-07 0.045 0.1128 0.9129 0.9000   1.00  1  1  3.4E-07
 11 :  -3.18E+00 4.48E-08 0.000 0.1329 0.9071 0.9000   1.00  2  2  5.3E-08
 12 :  -3.18E+00 4.46E-09 0.000 0.0994 0.9000 0.0000   1.00  2  2  1.2E-08

iter seconds digits       c*x               b*y
 12      0.1   Inf -3.1819805111e+00 -3.1819804667e+00
|Ax-b| =   1.1e-08, [Ay-c]_+ =   2.7E-09, |x|=  1.7e+01, |y|=  4.3e+00

Detailed timing (sec)
   Pre          IPM          Post
1.000E-02    5.000E-02    0.000E+00    
Max-norms: ||b||=1, ||c|| = 1,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 4319.78.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +3.18198
Nonzero lambda values for D design:
lambda(1) = 0.50001
lambda(10) = 0.49999
 
Calling sedumi: 32 variables, 11 equality constraints
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 11, order n = 27, dim = 39, blocks = 3
nnz(A) = 146 + 0, nnz(ADA) = 81, nnz(L) = 46
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            2.88E+02 0.000
  1 :   5.52E-01 7.08E+01 0.000 0.2457 0.9000 0.9000   1.01  1  1  3.7E+01
  2 :   8.89E-01 2.21E+01 0.000 0.3121 0.9000 0.9000   1.44  1  1  9.3E+00
  3 :   9.34E-01 7.50E+00 0.000 0.3392 0.9000 0.9000   1.47  1  1  2.5E+00
  4 :   9.08E-01 2.75E+00 0.000 0.3662 0.9000 0.9000   1.33  1  1  8.2E-01
  5 :   8.68E-01 7.45E-01 0.000 0.2715 0.9000 0.9000   1.04  1  1  2.2E-01
  6 :   8.68E-01 1.50E-02 0.319 0.0201 0.9900 0.0000   0.97  1  1  2.0E-02
  7 :   8.51E-01 7.40E-04 0.000 0.0493 0.9900 0.9900   0.99  1  1  9.8E-04
  8 :   8.50E-01 6.87E-05 0.066 0.0928 0.9069 0.9000   0.99  1  1  1.3E-04
  9 :   8.50E-01 3.63E-06 0.173 0.0529 0.9278 0.9000   0.99  1  1  1.6E-05
 10 :   8.50E-01 6.37E-07 0.099 0.1754 0.9194 0.9000   1.00  1  1  3.4E-06
 11 :   8.50E-01 1.32E-07 0.000 0.2081 0.9057 0.9000   1.00  1  1  7.3E-07
 12 :   8.50E-01 2.56E-08 0.000 0.1931 0.7982 0.9000   1.00  1  1  1.4E-07
 13 :   8.50E-01 2.02E-09 0.030 0.0791 0.9217 0.9900   1.00  1  1  1.3E-08

iter seconds digits       c*x               b*y
 13      0.2   Inf  8.4952798035e-01  8.4952803828e-01
|Ax-b| =   6.2e-09, [Ay-c]_+ =   8.9E-09, |x|=  8.1e+00, |y|=  1.7e+00

Detailed timing (sec)
   Pre          IPM          Post
1.000E-02    2.200E-01    0.000E+00    
Max-norms: ||b||=1, ||c|| = 1,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 19.0554.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.849528
Nonzero lambda values for A design:
lambda(1) = 0.29654
lambda(10) = 0.37799
lambda(20) = 0.32548
 
Calling sedumi: 23 variables, 3 equality constraints
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 3, order n = 23, dim = 25, blocks = 2
nnz(A) = 62 + 0, nnz(ADA) = 9, nnz(L) = 6
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            2.91E+01 0.000
  1 :  -3.67E+00 1.07E+01 0.000 0.3665 0.9000 0.9000   1.92  1  1  2.5E+01
  2 :  -1.62E+00 4.94E+00 0.000 0.4633 0.9000 0.9000   3.76  1  1  5.2E+00
  3 :  -1.78E+00 1.06E+00 0.000 0.2153 0.9000 0.9000   1.49  1  1  8.6E-01
  4 :  -1.80E+00 6.76E-02 0.000 0.0636 0.9900 0.9900   1.39  1  1  4.5E-02
  5 :  -1.80E+00 3.91E-05 0.000 0.0006 0.9999 0.9999   1.04  1  1  2.6E-05
  6 :  -1.80E+00 1.54E-12 0.487 0.0000 1.0000 1.0000   1.00  1  1  3.4E-12

iter seconds digits       c*x               b*y
  6      0.1   Inf -1.8000000000e+00 -1.8000000000e+00
|Ax-b| =   3.0e-12, [Ay-c]_+ =   0.0E+00, |x|=  8.2e-01, |y|=  2.1e+00

Detailed timing (sec)
   Pre          IPM          Post
2.000E-02    1.000E-01    1.000E-02    
Max-norms: ||b||=1, ||c|| = 9,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.09601.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1.8
Nonzero lambda values for E design:
lambda(10) = 0.2
lambda(20) = 0.8