math_utils.py

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