resources/libraries/python/PLRsearch/log_plus.py

   1 # Copyright (c) 2018 Cisco and/or its affiliates.
   2 # Licensed under the Apache License, Version 2.0 (the "License");
   3 # you may not use this file except in compliance with the License.
   4 # You may obtain a copy of the License at:
   5 #
   6 #     http://www.apache.org/licenses/LICENSE-2.0
   7 #
   8 # Unless required by applicable law or agreed to in writing, software
   9 # distributed under the License is distributed on an "AS IS" BASIS,
  10 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  11 # See the License for the specific language governing permissions and
  12 # limitations under the License.
  13
  14 """Module holding functions for avoiding rounding underflows.
  15
  16 Some applications wish to manipulate non-negative numbers
  17 which are many orders of magnitude apart.
  18 In those circumstances, it is useful to store
  19 the logarithm of the intended number.
  20
  21 As math.log(0.0) raises an exception (instead returning -inf or nan),
  22 and 0.0 might be a result of computation caused only by rounding error,
  23 functions of this module use None as -inf.
  24
  25 TODO: Figure out a more performant way of handling -inf.
  26
  27 The functions handle the common task of adding or substracting
  28 two numbers where both operands and the result is given in logarithm form.
  29 There are conditionals to make sure overflow does not happen (if possible)
  30 during the computation."""
  31
  32 import math
  33
  34
  35 def log_plus(first, second):
  36     """Return logarithm of the sum of two exponentials.
  37
  38     Basically math.log(math.exp(first) + math.exp(second))
  39     which avoids overflow and uses None as math.log(0.0).
  40
  41     TODO: replace with scipy.special.logsumexp ? Test it.
  42
  43     :param first: Logarithm of the first number to add (or None if zero).
  44     :param second: Logarithm of the second number to add (or None if zero).
  45     :type first: float
  46     :type second: float
  47     :returns: Logarithm of the sum (or None if zero).
  48     :rtype: float
  49     """
  50
  51     if first is None:
  52         return second
  53     if second is None:
  54         return first
  55     if second > first:
  56         return second + math.log(1.0 + math.exp(first - second))
  57     else:
  58         return first + math.log(1.0 + math.exp(second - first))
  59
  60
  61 def log_minus(first, second):
  62     """Return logarithm of the difference of two exponentials.
  63
  64     Basically math.log(math.exp(first) - math.exp(second))
  65     which avoids overflow and uses None as math.log(0.0).
  66
  67     TODO: Support zero difference?
  68     TODO: replace with scipy.special.logsumexp ? Test it.
  69
  70     :param first: Logarithm of the number to subtract from (or None if zero).
  71     :param second: Logarithm of the number to subtract (or None if zero).
  72     :type first: float
  73     :type second: float
  74     :returns: Logarithm of the difference.
  75     :rtype: float
  76     :raises RuntimeError: If the difference would be non-positive.
  77     """
  78
  79     if first is None:
  80         raise RuntimeError("log_minus: does not suport None first")
  81     if second is None:
  82         return first
  83     if second >= first:
  84         raise RuntimeError("log_minus: first has to be bigger than second")
  85     factor = -math.expm1(second - first)
  86     if factor <= 0.0:
  87         raise RuntimeError("log_minus: non-positive number to log")
  88     else:
  89         return first + math.log(factor)