# Source code for stonesoup.types.numeric

```# -*- coding: utf-8 -*-
from math import log, log1p, exp, trunc, ceil, floor
from numbers import Real, Integral

import numpy as np

[docs]class Probability(Real):
"""Probability class.

Similar to a float, but value stored as natural log value internally.
All operations are attempted with log values where possible, and failing
that a float will be returned instead.
"""

def __init__(self, value, *, log_value=False):
if log_value:
self._log_value = value
else:
try:
self._log_value = self._log(value)
except ValueError:
raise ValueError("value must be greater than 0")

@property
def log_value(self):
return self._log_value

@staticmethod
def _log(other):
if isinstance(other, Probability):
return other.log_value
elif other == 0:
return float("-inf")
else:
return log(other)

def __hash__(self):
value = float(self)
if value == 0 and self.log_value != float("-inf"):  # Too close to zero
# Add string so doesn't have same hash as the log value itself.
return hash(('log', self._log_value))
else:
return hash(value)

def __eq__(self, other):
if other < 0:
return False
return self.log_value == self._log(other)

def __le__(self, other):
if other < 0:
return False
return self.log_value <= self._log(other)

def __lt__(self, other):
if other < 0:
return False
return self.log_value < self._log(other)

def __ge__(self, other):
if other < 0:
return True
return self.log_value >= self._log(other)

def __gt__(self, other):
if other < 0:
return True
return self.log_value > self._log(other)

if other < 0:
return self - -other

log_other = self._log(other)
if self.log_value > log_other:
log_l, log_s = self.log_value, log_other
elif self.log_value < log_other:
log_l, log_s = log_other, self.log_value
else:  # Must be equal, so just double value
return self * 2

if log_s == float("-inf"):  # Just return largest value
return Probability(log_l, log_value=True)

return Probability(log_l + log1p(exp(log_s - log_l)),
log_value=True)

return self + other

def __sub__(self, other):
if other < 0:
return self + -other

log_other = self._log(other)
if self.log_value > log_other:
log_l, log_s = self.log_value, log_other
elif self.log_value < log_other:  # Result will be negative
return float(self) - other
else:  # Must be equal, so return 0
return Probability(float("-inf"), log_value=True)

if log_s == float("-inf"):  # Just return largest value
return Probability(log_l, log_value=True)

exp_diff = exp(log_s - log_l)
if exp_diff == 1:  # Diff too small, so result is effectively zero
return Probability(float("-inf"), log_value=True)

return Probability(log_l + log1p(-exp_diff),
log_value=True)

def __rsub__(self, other):
if other < 0:  # Result will be negative
return other + -float(self)

log_other = self._log(other)
if log_other > self.log_value:
log_l, log_s = log_other, self.log_value
elif log_other < self.log_value:  # Result will be negative
return other + -float(self)
else:  # Must be equal, so return 0
return Probability(float("-inf"), log_value=True)

if log_s == float("-inf"):  # Just return largest value
return Probability(log_l, log_value=True)

exp_diff = exp(log_s - log_l)
if exp_diff == 1:  # Diff too small, so result is effectively zero
return Probability(float("-inf"), log_value=True)

return Probability(log_l + log1p(-exp_diff),
log_value=True)

def __mul__(self, other):
try:
return Probability(self.log_value + self._log(other),
log_value=True)
except ValueError:
return float(self) * other

def __rmul__(self, other):
return self * other

def __truediv__(self, other):
try:
return Probability(self.log_value - self._log(other),
log_value=True)
except ValueError:
return float(self) / other

def __rtruediv__(self, other):
try:
return Probability(self._log(other) - self.log_value,
log_value=True)
except ValueError:
return other / float(self)

def __floordiv__(self, other):
return floor(self / other)

def __rfloordiv__(self, other):
return floor(other / self)

def __mod__(self, other):
try:
return Probability(float(self) % other)
except ValueError:
return float(self) % other

def __rmod__(self, other):
return Probability(other % float(self))

def __pow__(self, exponent):
if isinstance(exponent, Probability):
exponent = float(exponent)
return Probability(exponent * self.log_value, log_value=True)

def __rpow__(self, base):
return Probability(base ** float(self))

def __neg__(self):
return -float(self)

def __pos__(self):
return Probability(self)

def __abs__(self):
return Probability(self)

def __float__(self):
return exp(self.log_value)

def __round__(self, ndigits=None):
value = round(float(self), ndigits)
if isinstance(value, Integral):
return value
else:
return Probability(value)

def __trunc__(self):
return trunc(float(self))

def __floor__(self):
return floor(float(self))

def __ceil__(self):
return ceil(float(self))

def __repr__(self):
value = float(self)
if value == 0 and self.log_value != float("-inf"):  # Too close to zero
return "Probability({!r}, log_value=True)".format(self.log_value)
else:
return "Probability({!r})".format(float(self))

def __str__(self):
value = float(self)
if value == 0 and self.log_value != float("-inf"):  # Too close to zero
return "exp({})".format(self.log_value)
else:
return str(value)

[docs]    def sqrt(self):
"""Square root which can be called by NumPy"""
return self ** 0.5

[docs]    def log(self):
"""Log which can be called by NumPy"""
return self.log_value

[docs]    @classmethod
def sum(cls, values):
"""Carry out LogSumExp"""
log_values = np.array([cls._log(value) for value in values])
if len(log_values) == 0:
return Probability(0)

max_log_value = max(log_values)
# Check if all values are zero
if max_log_value == float('-inf'):
return Probability(0)

value_sum = np.sum(np.exp(log_values - max_log_value))

return Probability(cls._log(value_sum) + max_log_value, log_value=True)
```