math_utils.py

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