"""maxentutils.py: Utility routines for the maximum entropy module. Most
of them are either Python replacements for the corresponding Fortran
routines or wrappers around matrices to allow the maxent module to
manipulate ndarrays, scipy sparse matrices, and PySparse matrices a
common interface.
Perhaps the logsumexp() function belongs under the utils/ branch where other
modules can access it more easily.
Copyright: Ed Schofield, 2003-2006
License: BSD-style (see LICENSE.txt in main source directory)
"""
# Future imports must come before any code in 2.5
from __future__ import division
__author__ = "Ed Schofield"
__version__ = '2.0'
import random
import math
import cmath
import numpy
from numpy import log,exp,asarray,ndarray
from scipy import sparse
def logsumexp(a):
"""Compute the log of the sum of exponentials log(e^{a_1}+...e^{a_n})
of the components of the array a, avoiding numerical overflow.
"""
a = asarray(a)
a_max = a.max()
return a_max + log((exp(a-a_max)).sum())
def _logsumexpcomplex(values):
"""A version of logsumexp that should work if the values passed are
complex-numbered, such as the output of robustarraylog(). So we
expect:
cmath.exp(logsumexpcomplex(robustarraylog(values))) ~= sum(values,axis=0)
except for a small rounding error in both real and imag components.
The output is complex. (To recover just the real component, use
A.real, where A is the complex return value.)
"""
if len(values) == 0:
return 0.0
iterator = iter(values)
# Get the first element
while True:
# Loop until we have a value greater than -inf
try:
b_i = iterator.next() + 0j
except StopIteration:
# empty
return float('-inf')
if b_i.real != float('-inf'):
break
# Now the rest
for a_i in iterator:
a_i += 0j
if b_i.real > a_i.real:
increment = robustlog(1.+cmath.exp(a_i - b_i))
# print "Increment is " + str(increment)
b_i = b_i + increment
else:
increment = robustlog(1.+cmath.exp(b_i - a_i))
# print "Increment is " + str(increment)
b_i = a_i + increment
return b_i
def logsumexp_naive(values):
"""For testing logsumexp(). Subject to numerical overflow for large
values (e.g. 720).
"""
s = 0.0
for x in values:
s += math.exp(x)
return math.log(s)
def robustlog(x):
"""Returns log(x) if x > 0, the complex log cmath.log(x) if x < 0,
or float('-inf') if x == 0.
"""
if x == 0.:
return float('-inf')
elif type(x) is complex or (type(x) is float and x < 0):
return cmath.log(x)
else:
return math.log(x)
def _robustarraylog(x):
""" An array version of robustlog. Operates on a real array x.
"""
arraylog = emptyarray(len(x), numpy.Complex64)
for i in range(len(x)):
xi = x[i]
if xi > 0:
arraylog[i] = math.log(xi)
elif xi == 0.:
arraylog[i] = float('-inf')
else:
arraylog[i] = cmath.log(xi)
return arraylog
#try:
# from logsumexp import logsumexp, logsumexpcomplex, robustarraylog
#except:
# print "Warning: could not load the fast FORTRAN library for logsumexp()."
# logsumexp = _logsumexp
# logsumexpcomplex = _logsumexpcomplex
# robustarraylog = _robustarraylog
# pass
def arrayexp(x):
"""Returns the elementwise antilog of the real array x. We try to
exponentiate with numpy.exp() and, if that fails, with python's
math.exp(). numpy.exp() is about 10 times faster but throws an
OverflowError exception for numerical underflow (e.g. exp(-800),
whereas python's math.exp() just returns zero, which is much more
helpful.
"""
try:
ex = numpy.exp(x)
except OverflowError:
print "Warning: OverflowError using numpy.exp(). Using slower Python"\
" routines instead!"
ex = numpy.empty(len(x), float)
for j in range(len(x)):
ex[j] = math.exp(x[j])
return ex
def arrayexpcomplex(x):
"""Returns the elementwise antilog of the vector x. We try to
exponentiate with numpy.exp() and, if that fails, with python's
math.exp(). numpy.exp() is about 10 times faster but throws an
OverflowError exception for numerical underflow (e.g. exp(-800),
whereas python's math.exp() just returns zero, which is much more
helpful.
"""
try:
ex = numpy.exp(x).real
except OverflowError:
ex = numpy.empty(len(x), float)
try:
for j in range(len(x)):
ex[j] = math.exp(x[j])
except TypeError:
# Perhaps x[j] is complex. If so, try using the complex
# exponential and returning the real part.
for j in range(len(x)):
ex[j] = cmath.exp(x[j]).real
return ex
def sample_wr(population, k):
"""Chooses k random elements (with replacement) from a population.
(From the Python Cookbook).
"""
n = len(population)
_random, _int = random.random, int # speed hack
return [population[_int(_random() * n)] for i in xrange(k)]
def densefeatures(f, x):
"""Returns a dense array of non-zero evaluations of the functions fi
in the list f at the point x.
"""
return numpy.array([fi(x) for fi in f])
def densefeaturematrix(f, sample):
"""Returns an (m x n) dense array of non-zero evaluations of the
scalar functions fi in the list f at the points x_1,...,x_n in the
list sample.
"""
# Was: return numpy.array([[fi(x) for fi in f] for x in sample])
m = len(f)
n = len(sample)
F = numpy.empty((m, n), float)
for i in xrange(m):
f_i = f[i]
for j in xrange(n):
x = sample[j]
F[i,j] = f_i(x)
#for j in xrange(n):
# x = sample[j]
# for i in xrange(m):
# F[j,i] = f[i](x)
return F
def sparsefeatures(f, x, format='csc_matrix'):
""" Returns an Mx1 sparse matrix of non-zero evaluations of the
scalar functions f_1,...,f_m in the list f at the point x.
If format='ll_mat', the PySparse module (or a symlink to it) must be
available in the Python site-packages/ directory. A trimmed-down
version, patched for NumPy compatibility, is available in the SciPy
sandbox/pysparse directory.
"""
m = len(f)
if format == 'll_mat':
import spmatrix
sparsef = spmatrix.ll_mat(m, 1)
elif format in ('dok_matrix', 'csc_matrix', 'csr_matrix'):
sparsef = sparse.dok_matrix((m, 1))
for i in xrange(m):
f_i_x = f[i](x)
if f_i_x != 0:
sparsef[i, 0] = f_i_x
if format == 'csc_matrix':
print "Converting to CSC matrix ..."
return sparsef.tocsc()
elif format == 'csr_matrix':
print "Converting to CSR matrix ..."
return sparsef.tocsr()
else:
return sparsef
def sparsefeaturematrix(f, sample, format='csc_matrix'):
"""Returns an (m x n) sparse matrix of non-zero evaluations of the scalar
or vector functions f_1,...,f_m in the list f at the points
x_1,...,x_n in the sequence 'sample'.
If format='ll_mat', the PySparse module (or a symlink to it) must be
available in the Python site-packages/ directory. A trimmed-down
version, patched for NumPy compatibility, is available in the SciPy
sandbox/pysparse directory.
"""
m = len(f)
n = len(sample)
if format == 'll_mat':
import spmatrix
sparseF = spmatrix.ll_mat(m, n)
elif format in ('dok_matrix', 'csc_matrix', 'csr_matrix'):
sparseF = sparse.dok_matrix((m, n))
else:
raise ValueError, "sparse matrix format not recognized"
for i in xrange(m):
f_i = f[i]
for j in xrange(n):
x = sample[j]
f_i_x = f_i(x)
if f_i_x != 0:
sparseF[i,j] = f_i_x
if format == 'csc_matrix':
return sparseF.tocsc()
elif format == 'csr_matrix':
return sparseF.tocsr()
else:
return sparseF
def dotprod(u,v):
"""This is a wrapper around general dense or sparse dot products.
It is not necessary except as a common interface for supporting
ndarray, scipy spmatrix, and PySparse arrays.
Returns the dot product of the (1 x m) sparse array u with the
(m x 1) (dense) numpy array v.
"""
#print "Taking the dot product u.v, where"
#print "u has shape " + str(u.shape)
#print "v = " + str(v)
try:
dotprod = numpy.array([0.0]) # a 1x1 array. Required by spmatrix.
u.matvec(v, dotprod)
return dotprod[0] # extract the scalar
except AttributeError:
# Assume u is a dense array.
return numpy.dot(u,v)
def innerprod(A,v):
"""This is a wrapper around general dense or sparse dot products.
It is not necessary except as a common interface for supporting
ndarray, scipy spmatrix, and PySparse arrays.
Returns the inner product of the (m x n) dense or sparse matrix A
with the n-element dense array v. This is a wrapper for A.dot(v) for
dense arrays and spmatrix objects, and for A.matvec(v, result) for
PySparse matrices.
"""
# We assume A is sparse.
(m, n) = A.shape
vshape = v.shape
try:
(p,) = vshape
except ValueError:
(p, q) = vshape
if n != p:
raise TypeError, "matrix dimensions are incompatible"
if isinstance(v, ndarray):
try:
# See if A is sparse
A.matvec
except AttributeError:
# It looks like A is dense
return numpy.dot(A, v)
else:
# Assume A is sparse
if sparse.isspmatrix(A):
innerprod = A.matvec(v) # This returns a float32 type. Why???
return innerprod
else:
# Assume PySparse format
innerprod = numpy.empty(m, float)
A.matvec(v, innerprod)
return innerprod
elif sparse.isspmatrix(v):
return A * v
else:
raise TypeError, "unsupported types for inner product"
def innerprodtranspose(A,v):
"""This is a wrapper around general dense or sparse dot products.
It is not necessary except as a common interface for supporting
ndarray, scipy spmatrix, and PySparse arrays.
Computes A^T V, where A is a dense or sparse matrix and V is a numpy
array. If A is sparse, V must be a rank-1 array, not a matrix. This
function is efficient for large matrices A. This is a wrapper for
u.T.dot(v) for dense arrays and spmatrix objects, and for
u.matvec_transp(v, result) for pysparse matrices.
"""
(m, n) = A.shape
#pdb.set_trace()
if hasattr(A, 'matvec_transp'):
# A looks like a PySparse matrix
if len(v.shape) == 1:
innerprod = numpy.empty(n, float)
A.matvec_transp(v, innerprod)
else:
raise TypeError, "innerprodtranspose(A,v) requires that v be " \
"a vector (rank-1 dense array) if A is sparse."
return innerprod
elif sparse.isspmatrix(A):
return A.rmatvec(v).transpose()
else:
# Assume A is dense
if isinstance(v, numpy.ndarray):
# v is also dense
if len(v.shape) == 1:
# We can't transpose a rank-1 matrix into a row vector, so
# we reshape it.
vm = v.shape[0]
vcolumn = numpy.reshape(v, (1, vm))
x = numpy.dot(vcolumn, A)
return numpy.reshape(x, (n,))
else:
#(vm, vn) = v.shape
# Assume vm == m
x = numpy.dot(numpy.transpose(v), A)
return numpy.transpose(x)
else:
raise TypeError, "unsupported types for inner product"
def rowmeans(A):
"""This is a wrapper for general dense or sparse dot products. It is
only necessary as a common interface for supporting ndarray,
scipy spmatrix, and PySparse arrays.
Returns a dense (m x 1) vector representing the mean of the rows of A,
which be an (m x n) sparse or dense matrix.
>>> a = numpy.array([[1,2],[3,4]], float)
>>> rowmeans(a)
array([ 1.5, 3.5])
"""
if type(A) is numpy.ndarray:
return A.mean(1)
else:
# Assume it's sparse
try:
n = A.shape[1]
except AttributeError:
raise TypeError, \
"rowmeans() only works with sparse and dense arrays"
rowsum = innerprod(A, numpy.ones(n, float))
return rowsum / float(n)
def columnmeans(A):
"""This is a wrapper for general dense or sparse dot products. It is
only necessary as a common interface for supporting ndarray,
scipy spmatrix, and PySparse arrays.
Returns a dense (1 x n) vector with the column averages of A, which can
be an (m x n) sparse or dense matrix.
>>> a = numpy.array([[1,2],[3,4]],'d')
>>> columnmeans(a)
array([ 2., 3.])
"""
if type(A) is numpy.ndarray:
return A.mean(0)
else:
# Assume it's sparse
try:
m = A.shape[0]
except AttributeError:
raise TypeError, \
"columnmeans() only works with sparse and dense arrays"
columnsum = innerprodtranspose(A, numpy.ones(m, float))
return columnsum / float(m)
def columnvariances(A):
"""This is a wrapper for general dense or sparse dot products. It
is not necessary except as a common interface for supporting ndarray,
scipy spmatrix, and PySparse arrays.
Returns a dense (1 x n) vector with unbiased estimators for the column
variances for each column of the (m x n) sparse or dense matrix A. (The
normalization is by (m - 1).)
>>> a = numpy.array([[1,2], [3,4]], 'd')
>>> columnvariances(a)
array([ 2., 2.])
"""
if type(A) is numpy.ndarray:
return numpy.std(A,0)**2
else:
try:
m = A.shape[0]
except AttributeError:
raise TypeError, \
"columnvariances() only works with sparse and dense arrays"
means = columnmeans(A)
return columnmeans((A-means)**2) * (m/(m-1.0))
def flatten(a):
"""Flattens the sparse matrix or dense array/matrix 'a' into a
1-dimensional array
"""
if sparse.isspmatrix(a):
return a.A.flatten()
else:
return numpy.asarray(a).flatten()
class DivergenceError(Exception):
"""Exception raised if the entropy dual has no finite minimum.
"""
def __init__(self, message):
self.message = message
Exception.__init__(self)
def __str__(self):
return repr(self.message)
def _test():
import doctest
doctest.testmod()
if __name__ == "__main__":
_test()
|