% Generates figure 11.7 of Boyd & Vandenberghe, Convex Optimization
%
% Barrier method for standard form LPs
% 
%     minimize   c'*x 
%     subject to A*x = b
%                x >= 0

function figure_11_7()

% Figure 11.7.  Gap versus Newton iterations for three problems.

disp('m=50, n=100.')
m = 50; 
n = 100;
A = randn(m,n);  
x0 = rand(n,1);  
b = A*x0; 
z = randn(m,1);  
c = A'*z + rand(n,1);
[x, iters, gaps] = stdlp(A,b,c,x0,100,1e-3,1e-5);
l = length(gaps);  iters1 = [];  gaps1 = [];
for i=1:l-1
    iters1 = [iters1,  iters(i)-1,  iters(i+1)-1];  
    gaps1 = [gaps1, gaps(i), gaps(i)]; 
end;
iters1 = [iters1, iters(l)-1] - iters1(1); 
gaps1 = [gaps1, gaps(l)]; 

disp('m=500, n=1000.')
m = 500; 
n = 1000;
A = randn(m,n); 
x0 = rand(n,1); 
b = A*x0; 
z = rand(m,1);  
c = A'*z + rand(n,1);
[x, iters, gaps] = stdlp(A,b,c,x0,100,1e-3,1e-5);
l = length(gaps);  iters2 = [];  gaps2 = [];
for i=1:l-1
   iters2 = [iters2, iters(i)-1, iters(i+1)-1];  
   gaps2 = [gaps2, gaps(i), gaps(i)]; 
end;
iters2 = [iters2, iters(l)-1] - iters2(1);  
gaps2 = [gaps2, gaps(l)]; 

disp('m=1000, n=2000.')
m = 1000; 
n = 2000;
A = randn(m,n); 
x0 = rand(n,1); 
b = A*x0; 
z = rand(m,1);  
c = A'*z + rand(n,1);
[x, iters, gaps] = stdlp(A,b,c,x0,100,1e-3,1e-5);
l = length(gaps);  iters3 = [];  gaps3 = [];
for i=1:l-1
   iters3 = [iters3, iters(i)-1, iters(i+1)-1];  
   gaps3 = [gaps3, gaps(i), gaps(i)]; 
end;
iters3 = [iters3, iters(l)-1] - iters3(1);  
gaps3 = [gaps3, gaps(l)]; 


semilogy(iters1,gaps1, iters2,gaps2, '--', iters3,gaps3, '-');
text(iters1(length(iters1)), gaps1(length(iters1)),'m=50');
text(iters2(length(iters2)), gaps2(length(iters2)),'m=100');
text(iters3(length(iters3)), gaps3(length(iters3)),'m=500');
xlabel('Newton iterations'); ylabel('Duality gap');


function [x, inniters, gaps] = stdlp(A, b, c, x0, mu, tol, reltol)

% [x, inniters, gaps] = stdlp(A, b, c, x0, mu, tol, reltol)
%
% (primal) minimize  c'*x        (dual)  maximize   b'*z 
%          subjec to A*x = b             subject to A'*z <= c
%                    x >= 0
%
% x0:  strictly feasible starting point, not necesarily on central path
% inniters:  array with number of Newton iterations per outer iteration
% gaps:  array with duality gaps at the end of each outer iteration

MAXITERS = 500;   
ALPHA = 0.01;
BETA = 0.5;
NTTOL = 1e-5;     % stop inner iteration if lambda^2/2 < NTTOL

[m,n] = size(A);
t = 1;
x = x0;
gaps = [];
inniters = [];
for k=1:MAXITERS
    val = t*c'*x - sum(log(x));
    g = t*c - 1./x;
    z =  -(A*spdiags(x.^2,0,n,n)*A') \ (A*(g.*(x.^2)) );
    dx = - (g+A'*z) .* (x.^2);
    fprime = (g+A'*z)'*dx;
    if ((-fprime/2) < NTTOL)
        z = -z/t;
        dualobj = b'*z;
        gap = x'*(c-A'*z);  
        inniters = [inniters, k];  
        gaps = [gaps,gap];
        if (gap < tol),  return;  end;
        t = min(t*mu, (n+1)/tol);  
    else
        s = 1;
        while (min(x+s*dx) < 0), s = BETA*s; end;
        while (s*(t*c+A'*z)'*dx - sum(log(1+s*dx./x)) > ...
            ALPHA*s*fprime), 
            s = BETA*s; 
        end;
        x = x+s*dx;
    end;
end;