A numeric type for Python that performs nonnegative float computations in log-space.

logtype.py

__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