1 # Copyright (c) 2018 Cisco and/or its affiliates.
3 # you may not use this file except in compliance with the License.
4 # You may obtain a copy of the License at:
5 #
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.
14 """Module holding functions for avoiding rounding underflows.
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.
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.
25 TODO: Figure out a more performant way of handling -inf.
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."""
32 import math
35 def log_plus(first, second):
36     """Return logarithm of the sum of two exponentials.
38     Basically math.log(math.exp(first) + math.exp(second))
39     which avoids overflow and uses None as math.log(0.0).
41     TODO: replace with scipy.special.logsumexp ? Test it.
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     """
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))
61 def log_minus(first, second):
62     """Return logarithm of the difference of two exponentials.
64     Basically math.log(math.exp(first) - math.exp(second))
65     which avoids overflow and uses None as math.log(0.0).
67     TODO: Support zero difference?
68     TODO: replace with scipy.special.logsumexp ? Test it.
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     """
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)