#!/usr/bin/python

# The analytic centering example at the end of chapter 4 (The LAPACK 
# interface).

from cvxopt import matrix, log, mul, div, blas, lapack, base
from math import sqrt

def acent(A,b):
    """  
    Computes analytic center of A*x <= b with A m by n of rank n. 
    We assume that b > 0 and the feasible set is bounded.
    """

    MAXITERS = 100
    ALPHA = 0.01
    BETA = 0.5
    TOL = 1e-8

    ntdecrs = []
    m, n = A.size
    x = matrix(0.0, (n,1))
    H = matrix(0.0, (n,n))

    for iter in xrange(MAXITERS):
        
        # Gradient is g = A^T * (1./(b-A*x)).
        d = (b-A*x)**-1
        g = A.T * d

        # Hessian is H = A^T * diag(1./(b-A*x))^2 * A.
        Asc = mul( d[:,n*[0]], A)
        blas.syrk(Asc, H, trans='T')

        # Newton step is v = H^-1 * g.
        v = -g
        lapack.posv(H, v)

        # Directional derivative and Newton decrement.
	lam = blas.dot(g, v)
        ntdecrs += [ sqrt(-lam) ]
        print "%2d.  Newton decr. = %3.3e" %(iter,ntdecrs[-1])
        if ntdecrs[-1] < TOL: return x, ntdecrs

        # Backtracking line search.
        y = mul(A*v, d)
	step = 1.0
        while 1-step*max(y) < 0: step *= BETA 
	while True:
            if -sum(log(1-step*y)) < ALPHA*step*lam: break
	    step *= BETA
        x += step*v


# Generate an analytic centering problem  
#
#    -b1 <=  Ar*x <= b2 
#
# with random mxn Ar and random b1, b2.

m, n  = 500, 500
Ar = base.normal(m,n);
A = matrix([Ar, -Ar])
b = base.uniform(2*m,1)

x, ntdecrs = acent(A, b)  
try: 
    import pylab
except ImportError: 
    pass
else:
    pylab.semilogy(range(len(ntdecrs)), ntdecrs, 'o', 
             range(len(ntdecrs)), ntdecrs, '-')
    pylab.xlabel('Iteration number')
    pylab.ylabel('Newton decrement')
    pylab.show()