#!/usr/bin/python

# Figure 6.2, page 297.
# Penalty approximation.
#
# The problem data are not the same as in the book figure.

import pylab, numpy
from cvxopt import lapack, solvers, matrix, spdiag, log, div, normal
from cvxopt.modeling import variable, op, max, sum 
solvers.options['show_progress'] = 0

m, n = 100, 30
A = normal(m,n)
b = normal(m,1)
b /= (1.1 * max(abs(b)))   # Make x = 0 feasible for log barrier.


# l1 approximation
#
# minimize || A*x + b ||_1

x = variable(n)
op(sum(abs(A*x+b))).solve()
x1 = x.value

pylab.figure(1, facecolor='w', figsize=(10,10))
pylab.subplot(411)
nbins = 100
bins = [-1.5 + 3.0/(nbins-1)*k for k in xrange(nbins)]
pylab.hist( A*x1+b , numpy.array(bins))
nopts = 200
xs = -1.5 + 3.0/(nopts-1) * matrix(range(nopts))
pylab.plot(xs, (35.0/1.5) * abs(xs), 'g-')
pylab.axis([-1.5, 1.5, 0, 40])
pylab.ylabel('l1')
pylab.title('Penalty function approximation (fig. 6.2)')



# l2 approximation
#
# minimize || A*x + b ||_2

x = matrix(0.0, (m,1))
lapack.gels(+A, x)
x2 = x[:n]

pylab.subplot(412)
pylab.hist(A*x2+b, numpy.array(bins))
pylab.plot(xs, (8.0/1.5**2) * xs**2 , 'g-')
pylab.ylabel('l2')
pylab.axis([-1.5, 1.5, 0, 10])


# Deadzone approximation
#
# minimize sum(max(abs(A*x+b)-0.5, 0.0))

x = variable(n)
op(sum(max(abs(A*x+b)-0.5, 0.0))).solve()
xdz = x.value

pylab.subplot(413)
pylab.hist(A*xdz+b, numpy.array(bins))
pylab.plot(xs, 15.0/1.0 * matrix([ max(abs(xk)-0.5, 0.0) for xk 
    in xs ]), 'g-')
pylab.ylabel('Deadzone')
pylab.axis([-1.5, 1.5, 0, 20])


# Log barrier
#
# minimize -sum (log ( 1.0 - A*x+b)**2)

def F(x=None, z=None):
    if x is None: return 0, matrix(0.0, (n,1))
    y = A*x+b
    if max(abs(y)) >= 1.0: return None
    f = -sum(log(1.0 - y**2))
    gradf = 2.0 * A.T * div(y, 1-y**2)
    if z is None: return f, gradf.T
    H = A.T * spdiag(2.0 * z[0] * div( 1.0+y**2, (1.0 - y**2)**2 )) * A
    return f, gradf.T, H
xlb = solvers.cp(F)['x']

pylab.subplot(414)
pylab.hist(A*xlb+b, numpy.array(bins))
nopts = 200
pylab.plot(xs, (8.0/1.5**2) * xs**2, 'g--')
xs2 = -0.99999 + (2*0.99999 /(nopts-1)) * matrix(range(nopts))
pylab.plot(xs2, -3.0 * log(1.0 - abs(xs2)**2), 'g-')
pylab.ylabel('Log barrier')
pylab.xlabel('residual')
pylab.axis([-1.5, 1.5, 0, 10])
pylab.show()