"""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()