Skip to content

Instantly share code, notes, and snippets.

@bshillingford
Last active December 18, 2015 07:39
Show Gist options
  • Save bshillingford/5748168 to your computer and use it in GitHub Desktop.
Save bshillingford/5748168 to your computer and use it in GitHub Desktop.
A numeric type for Python that performs nonnegative float computations in log-space.
__author__ = 'brendan'
__all__ = ['logfloat', 'LogFloat']
import numpy as np # numpy's logaddexp correctly handles LOGZERO
from math import log as _log, exp as _exp
@np.vectorize
def logfloat(x):
"""
Constructs a LogFloat from x, or vectorized (via np.vectorize) for each element of x.
See also: `LogFloat`
"""
return LogFloat(x)
class LogFloat(object):
"""
Overloads operators for easy log-space computations.
Note that complex numbers are NOT used, so negative numbers are unsupported. In other words, this class is
optimized for probability calculations. Special symbol LOGZERO is used for log(0), as in:
http://bozeman.genome.washington.edu/compbio/mbt599_2006/hmm_scaling_revised.pdf .
"""
LOGZERO = float("-inf")
def __init__(self, normal=None, log_space=None):
"""Turns a normal-space float into a log-space float."""
if normal is not None:
if log_space is not None:
raise ValueError("do not specify both")
self.l = self._safe_log(normal)
elif log_space is not None:
self.l = log_space
else:
self.l = None
@classmethod
def from_log(cls, log_space):
return cls(log_space=log_space)
@staticmethod
def _safe_log(x):
if x < 0:
raise ValueError("log of negative")
elif x == 0:
return LogFloat.LOGZERO
else:
return _log(x)
def __repr__(self):
return "LogFloat(log_space = {:+f}, normal = {:e})".format(self.l, float(self))
def __str__(self):
return repr(self)
def to_float(self):
"""Returns the normal-space equivalent value of this."""
return _exp(self.l) # correctly handles LOGZERO
def __float__(self):
return self.to_float()
def __add__(self, other):
assert isinstance(other, LogFloat)
return self.from_log(np.logaddexp(self.l, other.l))
def __div__(self, other):
assert isinstance(other, LogFloat)
return self.from_log(self.l - other.l)
def __mul__(self, other):
assert isinstance(other, LogFloat)
return self.from_log(self.l + other.l)
def __pow__(self, power, modulo=None):
if modulo is not None:
raise TypeError("modulo argument not allowed unless all arguments are integers")
return self.from_log(self.l * power)
def __iadd__(self, other):
assert isinstance(other, LogFloat)
self.l = np.logaddexp(self.l, other.l)
def __idiv__(self, other):
assert isinstance(other, LogFloat)
self.l -= other.l
def __imul__(self, other):
assert isinstance(other, LogFloat)
self.l += other.l
def __ipow__(self, power, modulo=None):
if modulo is not None:
raise TypeError("modulo argument not allowed unless all arguments are integers")
self.l *= power
def __lt__(self, other):
assert isinstance(other, LogFloat)
return self.l < other.l
def __le__(self, other):
assert isinstance(other, LogFloat)
return self.l <= other.l
def __eq__(self, other):
assert isinstance(other, LogFloat)
return self.l == other.l
def __ne__(self, other):
assert isinstance(other, LogFloat)
return self.l != other.l
def __gt__(self, other):
assert isinstance(other, LogFloat)
return self.l > other.l
def __ge__(self, other):
assert isinstance(other, LogFloat)
return self.l >= other.l
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment