math_utils.py

   1 #!/usr/bin/env python3
   2
   3 """Mathematical helpers."""
   4
   5 import functools
   6 import math
   7 from heapq import heappop, heappush
   8 from typing import List, Optional
   9
  10
  11 class NumericPopulation(object):
  12     """A numeric population with some statistics such as median, mean, pN,
  13     stdev, etc...
  14
  15     >>> pop = NumericPopulation()
  16     >>> pop.add_number(1)
  17     >>> pop.add_number(10)
  18     >>> pop.add_number(3)
  19     >>> pop.get_median()
  20     3
  21     >>> pop.add_number(7)
  22     >>> pop.add_number(5)
  23     >>> pop.get_median()
  24     5
  25     >>> pop.get_mean()
  26     5.2
  27     >>> round(pop.get_stdev(), 2)
  28     6.99
  29     >>> pop.get_percentile(20)
  30     3
  31     >>> pop.get_percentile(60)
  32     7
  33
  34     """
  35
  36     def __init__(self):
  37         self.lowers, self.highers = [], []
  38         self.aggregate = 0.0
  39         self.sorted_copy: Optional[List[float]] = None
  40
  41     def add_number(self, number: float):
  42         """O(2 log2 n)"""
  43
  44         if not self.highers or number > self.highers[0]:
  45             heappush(self.highers, number)
  46         else:
  47             heappush(self.lowers, -number)  # for lowers we need a max heap
  48         self.aggregate += number
  49         self._rebalance()
  50
  51     def _rebalance(self):
  52         if len(self.lowers) - len(self.highers) > 1:
  53             heappush(self.highers, -heappop(self.lowers))
  54         elif len(self.highers) - len(self.lowers) > 1:
  55             heappush(self.lowers, -heappop(self.highers))
  56
  57     def get_median(self) -> float:
  58         """Returns the approximate median (p50) so far in O(1) time."""
  59
  60         if len(self.lowers) == len(self.highers):
  61             return -self.lowers[0]
  62         elif len(self.lowers) > len(self.highers):
  63             return -self.lowers[0]
  64         else:
  65             return self.highers[0]
  66
  67     def get_mean(self) -> float:
  68         """Returns the mean (arithmetic mean) so far in O(1) time."""
  69
  70         count = len(self.lowers) + len(self.highers)
  71         return self.aggregate / count
  72
  73     def get_stdev(self) -> float:
  74         """Returns the stdev so far in O(n) time."""
  75
  76         mean = self.get_mean()
  77         variance = 0.0
  78         for n in self.lowers:
  79             n = -n
  80             variance += (n - mean) ** 2
  81         for n in self.highers:
  82             variance += (n - mean) ** 2
  83         return math.sqrt(variance)
  84
  85     def get_percentile(self, n: float) -> float:
  86         """Returns the number at approximately pn% (i.e. the nth percentile)
  87         of the distribution in O(n log n) time (expensive, requires a
  88         complete sort).  Not thread safe.  Caching does across
  89         multiple calls without an invocation to add_number.
  90
  91         """
  92         if n == 50:
  93             return self.get_median()
  94         count = len(self.lowers) + len(self.highers)
  95         if self.sorted_copy is not None:
  96             if count == len(self.sorted_copy):
  97                 index = round(count * (n / 100.0))
  98                 assert 0 <= index < count
  99                 return self.sorted_copy[index]
 100         self.sorted_copy = [-x for x in self.lowers]
 101         for x in self.highers:
 102             self.sorted_copy.append(x)
 103         self.sorted_copy = sorted(self.sorted_copy)
 104         index = round(count * (n / 100.0))
 105         assert 0 <= index < count
 106         return self.sorted_copy[index]
 107
 108
 109 def gcd_floats(a: float, b: float) -> float:
 110     if a < b:
 111         return gcd_floats(b, a)
 112
 113     # base case
 114     if abs(b) < 0.001:
 115         return a
 116     return gcd_floats(b, a - math.floor(a / b) * b)
 117
 118
 119 def gcd_float_sequence(lst: List[float]) -> float:
 120     if len(lst) <= 0:
 121         raise ValueError("Need at least one number")
 122     elif len(lst) == 1:
 123         return lst[0]
 124     assert len(lst) >= 2
 125     gcd = gcd_floats(lst[0], lst[1])
 126     for i in range(2, len(lst)):
 127         gcd = gcd_floats(gcd, lst[i])
 128     return gcd
 129
 130
 131 def truncate_float(n: float, decimals: int = 2):
 132     """
 133     Truncate a float to a particular number of decimals.
 134
 135     >>> truncate_float(3.1415927, 3)
 136     3.141
 137
 138     """
 139     assert 0 < decimals < 10
 140     multiplier = 10**decimals
 141     return int(n * multiplier) / multiplier
 142
 143
 144 def percentage_to_multiplier(percent: float) -> float:
 145     """Given a percentage (e.g. 155%), return a factor needed to scale a
 146     number by that percentage.
 147
 148     >>> percentage_to_multiplier(155)
 149     2.55
 150     >>> percentage_to_multiplier(45)
 151     1.45
 152     >>> percentage_to_multiplier(-25)
 153     0.75
 154
 155     """
 156     multiplier = percent / 100
 157     multiplier += 1.0
 158     return multiplier
 159
 160
 161 def multiplier_to_percent(multiplier: float) -> float:
 162     """Convert a multiplicative factor into a percent change.
 163
 164     >>> multiplier_to_percent(0.75)
 165     -25.0
 166     >>> multiplier_to_percent(1.0)
 167     0.0
 168     >>> multiplier_to_percent(1.99)
 169     99.0
 170
 171     """
 172     percent = multiplier
 173     if percent > 0.0:
 174         percent -= 1.0
 175     else:
 176         percent = 1.0 - percent
 177     percent *= 100.0
 178     return percent
 179
 180
 181 @functools.lru_cache(maxsize=1024, typed=True)
 182 def is_prime(n: int) -> bool:
 183     """
 184     Returns True if n is prime and False otherwise.  Obviously(?) very slow for
 185     very large input numbers.
 186
 187     >>> is_prime(13)
 188     True
 189     >>> is_prime(22)
 190     False
 191     >>> is_prime(51602981)
 192     True
 193
 194     """
 195     if not isinstance(n, int):
 196         raise TypeError("argument passed to is_prime is not of 'int' type")
 197
 198     # Corner cases
 199     if n <= 1:
 200         return False
 201     if n <= 3:
 202         return True
 203
 204     # This is checked so that we can skip middle five numbers in below
 205     # loop
 206     if n % 2 == 0 or n % 3 == 0:
 207         return False
 208
 209     i = 5
 210     while i * i <= n:
 211         if n % i == 0 or n % (i + 2) == 0:
 212             return False
 213         i = i + 6
 214     return True
 215
 216
 217 if __name__ == '__main__':
 218     import doctest
 219
 220     doctest.testmod()