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     1.75
  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         count = len(self.lowers) + len(self.highers) - 1
  84         return math.sqrt(variance) / count
  85
  86     def get_percentile(self, n: float) -> float:
  87         """Returns the number at approximately pn% (i.e. the nth percentile)
  88         of the distribution in O(n log n) time (expensive, requires a
  89         complete sort).  Not thread safe.  Caching does across
  90         multiple calls without an invocation to add_number.
  91
  92         """
  93         if n == 50:
  94             return self.get_median()
  95         count = len(self.lowers) + len(self.highers)
  96         if self.sorted_copy is not None:
  97             if count == len(self.sorted_copy):
  98                 index = round(count * (n / 100.0))
  99                 assert 0 <= index < count
 100                 return self.sorted_copy[index]
 101         self.sorted_copy = [-x for x in self.lowers]
 102         for x in self.highers:
 103             self.sorted_copy.append(x)
 104         self.sorted_copy = sorted(self.sorted_copy)
 105         index = round(count * (n / 100.0))
 106         assert 0 <= index < count
 107         return self.sorted_copy[index]
 108
 109
 110 def gcd_floats(a: float, b: float) -> float:
 111     if a < b:
 112         return gcd_floats(b, a)
 113
 114     # base case
 115     if abs(b) < 0.001:
 116         return a
 117     return gcd_floats(b, a - math.floor(a / b) * b)
 118
 119
 120 def gcd_float_sequence(lst: List[float]) -> float:
 121     if len(lst) <= 0:
 122         raise ValueError("Need at least one number")
 123     elif len(lst) == 1:
 124         return lst[0]
 125     assert len(lst) >= 2
 126     gcd = gcd_floats(lst[0], lst[1])
 127     for i in range(2, len(lst)):
 128         gcd = gcd_floats(gcd, lst[i])
 129     return gcd
 130
 131
 132 def truncate_float(n: float, decimals: int = 2):
 133     """
 134     Truncate a float to a particular number of decimals.
 135
 136     >>> truncate_float(3.1415927, 3)
 137     3.141
 138
 139     """
 140     assert 0 < decimals < 10
 141     multiplier = 10**decimals
 142     return int(n * multiplier) / multiplier
 143
 144
 145 def percentage_to_multiplier(percent: float) -> float:
 146     """Given a percentage (e.g. 155%), return a factor needed to scale a
 147     number by that percentage.
 148
 149     >>> percentage_to_multiplier(155)
 150     2.55
 151     >>> percentage_to_multiplier(45)
 152     1.45
 153     >>> percentage_to_multiplier(-25)
 154     0.75
 155
 156     """
 157     multiplier = percent / 100
 158     multiplier += 1.0
 159     return multiplier
 160
 161
 162 def multiplier_to_percent(multiplier: float) -> float:
 163     """Convert a multiplicative factor into a percent change.
 164
 165     >>> multiplier_to_percent(0.75)
 166     -25.0
 167     >>> multiplier_to_percent(1.0)
 168     0.0
 169     >>> multiplier_to_percent(1.99)
 170     99.0
 171
 172     """
 173     percent = multiplier
 174     if percent > 0.0:
 175         percent -= 1.0
 176     else:
 177         percent = 1.0 - percent
 178     percent *= 100.0
 179     return percent
 180
 181
 182 @functools.lru_cache(maxsize=1024, typed=True)
 183 def is_prime(n: int) -> bool:
 184     """
 185     Returns True if n is prime and False otherwise.  Obviously(?) very slow for
 186     very large input numbers.
 187
 188     >>> is_prime(13)
 189     True
 190     >>> is_prime(22)
 191     False
 192     >>> is_prime(51602981)
 193     True
 194
 195     """
 196     if not isinstance(n, int):
 197         raise TypeError("argument passed to is_prime is not of 'int' type")
 198
 199     # Corner cases
 200     if n <= 1:
 201         return False
 202     if n <= 3:
 203         return True
 204
 205     # This is checked so that we can skip middle five numbers in below
 206     # loop
 207     if n % 2 == 0 or n % 3 == 0:
 208         return False
 209
 210     i = 5
 211     while i * i <= n:
 212         if n % i == 0 or n % (i + 2) == 0:
 213             return False
 214         i = i + 6
 215     return True
 216
 217
 218 if __name__ == '__main__':
 219     import doctest
 220
 221     doctest.testmod()