#!/usr/bin/python

# Figure 7.8, page 384.
# Chernoff lower bound.

import pylab
from cvxopt import matrix, mul, exp, normal, solvers, blas
solvers.options['show_progress'] = False

# Extreme points and inequality description of Voronoi region around 
# first symbol (at the origin).
m = 6
V = matrix([ 1.0,  1.0, 
            -1.0,  2.0,
            -2.0,  1.0,
            -2.0, -1.0,
             0.0, -2.0,
             1.5, -1.0,
             1.0,  1.0 ], (2,m+1))

# A and b are lists with the inequality descriptions of the regions.
A = [ matrix( [-(V[1,:m] - V[1,1:]), V[0,:m] - V[0,1:]] ).T ]
b = [ mul(A[0], V[:,:m].T) * matrix(1.0, (2,1)) ]

# List of symbols.
C = [ matrix(0.0, (2,1)) ] + \
    [ 2.0 * b[0][k] / blas.nrm2(A[0][k,:])**2 * A[0][k,:].T for k in 
    xrange(m) ]

# Voronoi set around C[1]
A += [ matrix(0.0, (3,2)) ] 
b += [ matrix(0.0, (3,1)) ]
A[1][0,:] = -A[0][0,:]
b[1][0] = -b[0][0]
A[1][1,:] = (C[m] - C[1]).T
b[1][1] = 0.5 * A[1][1,:] * ( C[m] + C[1] )
A[1][2,:] = (C[2] - C[1]).T
b[1][2] = 0.5 * A[1][2,:] * ( C[2] + C[1] )

# Voronoi set around C[2], ..., C[5]
for k in xrange(2, 6):
    A += [ matrix(0.0, (3,2)) ] 
    b += [ matrix(0.0, (3,1)) ]
    A[k][0,:] = -A[0][k-1,:]
    b[k][0] = -b[0][k-1]
    A[k][1,:] = (C[k-1] - C[k]).T
    b[k][1] = 0.5 * A[k][1,:] * ( C[k-1] + C[k] )
    A[k][2,:] = (C[k+1] - C[k]).T
    b[k][2] = 0.5 * A[k][2,:] * ( C[k+1] + C[k] )

# Voronoi set around C[6]
A += [ matrix(0.0, (3,2)) ] 
b += [ matrix(0.0, (3,1)) ]
A[6][0,:] = -A[0][5,:]
b[6][0] = -b[0][5]
A[6][1,:] = (C[1] - C[6]).T
b[6][1] = 0.5 * A[6][1,:] * ( C[1] + C[6] )
A[6][2,:] = (C[5] - C[6]).T
b[6][2] = 0.5 * A[6][2,:] * ( C[5] + C[6] )


# For regions k=1, ..., 6, let pk be the optimal value of 
#
#     minimize    x'*x 
#     subject to  A*x <= b.
#
# The Chernoff upper bound is  1.0 - sum exp( - pk / (2 sigma^2)).

P = matrix([1.0, 0.0, 0.0, 1.0], (2,2))
q = matrix(0.0, (2,1))
optvals = matrix([ blas.nrm2( solvers.qp(P, q, A[k], b[k] )['x'] )**2
    for k in xrange(1,7) ])
nopts = 200
sigmas = 0.2 + (0.5 - 0.2)/nopts * matrix(range(nopts), tc='d')
bnds = [ 1.0 - sum( exp( - optvals / (2*sigma**2) )) for sigma 
    in sigmas]

pylab.figure(facecolor='w')
pylab.plot(sigmas, bnds, '-')
pylab.axis([0.2, 0.5, 0.9, 1.0])
pylab.title('Chernoff lower bound (fig. 7.8)')
pylab.xlabel('sigma')
pylab.ylabel('Probability of correct detection')
pylab.show()

if 0:  # uncomment out for the Monte Carlo estimation.
    N = 100000
    mcest = []
    ones = matrix(1.0, (1,m))
    for sigma in sigmas: 
        X = sigma * normal(2, N)
        S = b[0][:,N*[0]] - A[0]*X 
        S = ones * (S - abs(S))
        mcest += [ N - len(filter(lambda x: x < 0.0, S)) ]

    pylab.figure(facecolor='w')
    pylab.plot(sigmas, bnds, '-', sigmas, (1.0/N)*matrix(mcest), '--')
    pylab.plot(sigmas, bnds, '-')
    pylab.axis([0.2, 0.5, 0.9, 1.0])
    pylab.title('Chernoff lower bound (fig. 7.8)')
    pylab.xlabel('sigma')
    pylab.ylabel('Probability of correct detection')
    pylab.show()