diff --git a/README.md b/README.md index fd596dd..368d667 100644 --- a/README.md +++ b/README.md @@ -33,7 +33,7 @@ python -m doctest algorithm_lab/*.py Optional installation: `python -m pip install -e .`. For linting and formatting, install `.[dev]` and run `ruff check .` and `ruff format --check .`. -## Catalog (148 algorithms) +## Catalog (156 algorithms) | Module | Topic | What it demonstrates | | --- | --- | --- | @@ -185,6 +185,14 @@ install `.[dev]` and run `ruff check .` and `ruff format --check .`. | [interval_intersection](algorithm_lab/interval_intersection.py) | Intervals | Intersect two sorted disjoint lists of closed integer intervals | | [stable_matching](algorithm_lab/stable_matching.py) | Graphs | Gale-Shapley matching with complete strict preferences | | [zero_one_bfs](algorithm_lab/zero_one_bfs.py) | Graphs | Shortest distances in directed graphs with zero-or-one edge costs | +| [percentile](algorithm_lab/percentile.py) | Statistics | Linearly interpolated sample quantiles | +| [weighted_median](algorithm_lab/weighted_median.py) | Statistics | Lower weighted medians with exact integer weights | +| [entropy](algorithm_lab/entropy.py) | Statistics | Shannon entropy in bits from nonnegative weights | +| [jensen_shannon_divergence](algorithm_lab/jensen_shannon_divergence.py) | Statistics | Symmetric Jensen-Shannon divergence in bits | +| [softmax](algorithm_lab/softmax.py) | Statistics | Stable softmax by subtracting the maximum logit | +| [isotonic_regression](algorithm_lab/isotonic_regression.py) | Statistics | Exact unweighted isotonic regression with pooled adjacent violators | +| [kendall_tau](algorithm_lab/kendall_tau.py) | Statistics | Tie-adjusted Kendall rank correlation tau-b | +| [online_covariance](algorithm_lab/online_covariance.py) | Statistics | One-pass covariance for paired numeric observations | ## Validation diff --git a/algorithm_lab/_numeric.py b/algorithm_lab/_numeric.py new file mode 100644 index 0000000..25f6ef8 --- /dev/null +++ b/algorithm_lab/_numeric.py @@ -0,0 +1,25 @@ +"""Finite floating-point validation and scale-safe probability normalization.""" + +import math + + +def finite(value): + if isinstance(value, bool) or not isinstance(value, (int, float)): + raise ValueError("values must be finite integers or floats") + try: + result = float(value) + except OverflowError as exc: + raise ValueError("value exceeds the floating-point range") from exc + if not math.isfinite(result): + raise ValueError("values must be finite") + return result + + +def probabilities(weights): + values = [finite(value) for value in weights] + if not values or any(value < 0 for value in values) or max(values) == 0: + raise ValueError("weights must be nonnegative with positive total mass") + scale = max(values) + scaled = [value / scale for value in values] + total = math.fsum(scaled) + return [value / total for value in scaled] diff --git a/algorithm_lab/entropy.py b/algorithm_lab/entropy.py new file mode 100644 index 0000000..070f9f0 --- /dev/null +++ b/algorithm_lab/entropy.py @@ -0,0 +1,18 @@ +"""Shannon entropy in bits from nonnegative weights. + +Accept finite int/float weights, excluding bool, with positive total mass. +Normalize by the largest weight before summing to avoid overflow. Zero-mass +categories contribute zero. Return entropy in bits. O(n) time and space. +Extremely small probabilities may underflow to zero in floating point. + +>>> entropy([1, 1, 1, 1]) +2.0 +""" + +from math import fsum, log2 + +from algorithm_lab._numeric import probabilities + + +def entropy(weights): + return fsum(-p * log2(p) for p in probabilities(weights) if p > 0) diff --git a/algorithm_lab/isotonic_regression.py b/algorithm_lab/isotonic_regression.py new file mode 100644 index 0000000..31839d9 --- /dev/null +++ b/algorithm_lab/isotonic_regression.py @@ -0,0 +1,25 @@ +"""Exact unweighted isotonic regression with pooled adjacent violators. + +Accept int/Fraction values, excluding bool and float. Return nondecreasing +Fraction fitted values minimizing the sum of squared residuals. Adjacent +blocks merge only for a strict decrease. Empty input returns []. O(n) +rational operations and O(n) space; no observations are reordered. + +>>> isotonic_regression([3, 1, 2, 5]) +[Fraction(2, 1), Fraction(2, 1), Fraction(2, 1), Fraction(5, 1)] +""" + +from fractions import Fraction + + +def isotonic_regression(values): + blocks = [] + for value in values: + if isinstance(value, bool) or not isinstance(value, (int, Fraction)): + raise ValueError("values must be integers or Fractions") + blocks.append((Fraction(value), 1)) + while len(blocks) > 1 and blocks[-2][0] * blocks[-1][1] > blocks[-1][0] * blocks[-2][1]: + right, count_right = blocks.pop() + left, count_left = blocks.pop() + blocks.append((left + right, count_left + count_right)) + return [total / count for total, count in blocks for _ in range(count)] diff --git a/algorithm_lab/jensen_shannon_divergence.py b/algorithm_lab/jensen_shannon_divergence.py new file mode 100644 index 0000000..a512d14 --- /dev/null +++ b/algorithm_lab/jensen_shannon_divergence.py @@ -0,0 +1,28 @@ +"""Symmetric Jensen-Shannon divergence in bits. + +Inputs are equal-length vectors of finite nonnegative int/float weights, +excluding bool, each with positive mass. Normalize each independently and +compare against their equally weighted mixture. Return divergence in [0,1] +up to floating-point rounding; this is not the square-root distance. +O(n) time and space. Tiny probabilities may underflow to zero. + +>>> jensen_shannon_divergence([1, 0], [0, 1]) +1.0 +""" + +from math import fsum, log2 + +from algorithm_lab._numeric import probabilities + + +def jensen_shannon_divergence(left, right): + left, right = probabilities(left), probabilities(right) + if len(left) != len(right): + raise ValueError("distributions must have equal length") + terms = [] + for p, q in zip(left, right, strict=True): + if p: + terms.append(p * (log2(p) - log2(p + q) + 1) / 2) + if q: + terms.append(q * (log2(q) - log2(p + q) + 1) / 2) + return max(0.0, min(1.0, fsum(terms))) diff --git a/algorithm_lab/kendall_tau.py b/algorithm_lab/kendall_tau.py new file mode 100644 index 0000000..d6d7d21 --- /dev/null +++ b/algorithm_lab/kendall_tau.py @@ -0,0 +1,31 @@ +"""Tie-adjusted Kendall rank correlation tau-b. + +Inputs are equal-length finite int/float observations, excluding bool, with +at least two pairs. Return tau-b, accounting separately for ties on each +axis. A constant axis raises ValueError because the denominator is zero. +O(n^2) comparisons and O(n) copied storage. Ordering uses validated floats. + +>>> kendall_tau([1, 2, 3], [3, 2, 1]) +-1.0 +""" + +from math import sqrt + +from algorithm_lab._numeric import finite + + +def kendall_tau(left, right): + left, right = [finite(x) for x in left], [finite(x) for x in right] + if len(left) != len(right) or len(left) < 2: + raise ValueError("at least two paired observations are required") + concordance, untied_left, untied_right = 0, 0, 0 + for i in range(len(left)): + for j in range(i): + a = (left[i] > left[j]) - (left[i] < left[j]) + b = (right[i] > right[j]) - (right[i] < right[j]) + concordance += a * b + untied_left += a != 0 + untied_right += b != 0 + if not untied_left or not untied_right: + raise ValueError("correlation is undefined for a constant axis") + return concordance / sqrt(untied_left * untied_right) diff --git a/algorithm_lab/online_covariance.py b/algorithm_lab/online_covariance.py new file mode 100644 index 0000000..5251ecf --- /dev/null +++ b/algorithm_lab/online_covariance.py @@ -0,0 +1,34 @@ +"""One-pass covariance for paired numeric observations. + +observations yields (x,y) pairs of finite int/float values, excluding bool. +sample is a bool selecting n-1 rather than n as divisor. Require one pair +for population or two for sample covariance. Return a float; unrepresentable +intermediate state raises ValueError. O(n) time and O(1) state using a +paired Welford recurrence; observations need not fit in memory. + +>>> online_covariance([(1, 2), (2, 4), (3, 6)], sample=True) +2.0 +""" + +from math import isfinite + +from algorithm_lab._numeric import finite + + +def online_covariance(observations, *, sample=False): + if not isinstance(sample, bool): + raise ValueError("sample must be a boolean") + count, mean_x, mean_y, moment = 0, 0.0, 0.0, 0.0 + for x, y in observations: + x, y = finite(x), finite(y) + count += 1 + delta_x, delta_y = x - mean_x, y - mean_y + mean_x += delta_x / count + mean_y += delta_y / count + moment += delta_x * (y - mean_y) + if not all(isfinite(v) for v in (mean_x, mean_y, moment)): + raise ValueError("covariance exceeds floating-point range") + denominator = count - int(sample) + if denominator <= 0: + raise ValueError("not enough observations") + return moment / denominator diff --git a/algorithm_lab/percentile.py b/algorithm_lab/percentile.py new file mode 100644 index 0000000..b0a146c --- /dev/null +++ b/algorithm_lab/percentile.py @@ -0,0 +1,27 @@ +"""Linearly interpolated sample quantiles. + +Values are nonempty finite int/float observations excluding bool. q is a +finite number in [0,1], excluding bool. Return a float at position +q*(n-1) in sorted data (the common type-7 convention). O(n log n) time and +O(n) space. Input order is irrelevant and input is unchanged. + +>>> percentile([0, 10, 20, 30], 0.25) +7.5 +""" + +from math import floor + +from algorithm_lab._numeric import finite + + +def percentile(values, q): + values = sorted(finite(value) for value in values) + q = finite(q) + if not values or not 0 <= q <= 1: + raise ValueError("nonempty values and q in [0,1] are required") + position = q * (len(values) - 1) + index = floor(position) + if index == len(values) - 1: + return values[index] + fraction = position - index + return (1 - fraction) * values[index] + fraction * values[index + 1] diff --git a/algorithm_lab/softmax.py b/algorithm_lab/softmax.py new file mode 100644 index 0000000..7b6928f --- /dev/null +++ b/algorithm_lab/softmax.py @@ -0,0 +1,24 @@ +"""Stable softmax by subtracting the maximum logit. + +Accept a nonempty iterable of finite int/float logits, excluding bool. +Return nonnegative probabilities summing to approximately 1. Subtracting +the maximum avoids exponential overflow; very small tails can underflow +to zero. O(n) time and space, leaving input unchanged. + +>>> softmax([1000, 1000]) +[0.5, 0.5] +""" + +from math import exp, fsum + +from algorithm_lab._numeric import finite + + +def softmax(logits): + logits = [finite(value) for value in logits] + if not logits: + raise ValueError("logits must be nonempty") + maximum = max(logits) + values = [exp(value - maximum) for value in logits] + total = fsum(values) + return [value / total for value in values] diff --git a/algorithm_lab/weighted_median.py b/algorithm_lab/weighted_median.py new file mode 100644 index 0000000..bc293e4 --- /dev/null +++ b/algorithm_lab/weighted_median.py @@ -0,0 +1,31 @@ +"""Lower weighted medians with exact integer weights. + +Values and weights are equally sized integer iterables, excluding bool. +Weights are nonnegative with positive total; zero weights have no influence. +Return the smallest observed value whose cumulative weight reaches half the +total. This minimizes weighted absolute deviation. O(n log n) time, O(n) space. + +>>> weighted_median([10, 1, 5], [1, 2, 1]) +1 +""" + +from algorithm_lab._validation import integer + + +def weighted_median(values, weights): + values, weights = list(values), list(weights) + if len(values) != len(weights): + raise ValueError("values and weights must have equal length") + for value in values: + integer(value) + for weight in weights: + integer(weight, "weight", minimum=0) + total = sum(weights) + if not total: + raise ValueError("total weight must be positive") + cumulative = 0 + for value, weight in sorted(zip(values, weights, strict=True)): + cumulative += weight + if cumulative * 2 >= total: + return value + raise AssertionError("positive total must produce a median") diff --git a/tests/test_entropy.py b/tests/test_entropy.py new file mode 100644 index 0000000..d6b342b --- /dev/null +++ b/tests/test_entropy.py @@ -0,0 +1,19 @@ +import math +import unittest + +from algorithm_lab.entropy import entropy + + +class Tests(unittest.TestCase): + def test_known_distributions_and_scaling(self): + for n in range(1, 30): + self.assertAlmostEqual(entropy([1] * n), math.log2(n)) + self.assertEqual(entropy([0, 8, 0]), 0) + self.assertAlmostEqual(entropy([1, 3]), 0.8112781244591328) + self.assertAlmostEqual(entropy([1e308, 1e308]), 1) + self.assertAlmostEqual(entropy([1, 3]), entropy([1e100, 3e100])) + + def test_invalid(self): + for weights in [[], [0, 0], [-1, 2], [math.nan], [math.inf], [True], [10**1000]]: + with self.assertRaises(ValueError): + entropy(weights) diff --git a/tests/test_isotonic_regression.py b/tests/test_isotonic_regression.py new file mode 100644 index 0000000..c72275b --- /dev/null +++ b/tests/test_isotonic_regression.py @@ -0,0 +1,32 @@ +import unittest +from fractions import Fraction +from itertools import product + +from algorithm_lab.isotonic_regression import isotonic_regression + + +class Tests(unittest.TestCase): + def test_all_contiguous_block_partitions(self): + self.assertEqual(isotonic_regression([]), []) + for n in range(1, 6): + for values in product(range(3), repeat=n): + costs = [] + for mask in range(1 << (n - 1)): + ends = [0] + [i for i in range(1, n) if mask >> (i - 1) & 1] + [n] + candidate = [] + for a, b in zip(ends, ends[1:], strict=False): + candidate += [Fraction(sum(values[a:b]), b - a)] * (b - a) + if candidate == sorted(candidate): + costs.append( + sum((a - b) ** 2 for a, b in zip(values, candidate, strict=True)) + ) + actual = isotonic_regression(values) + self.assertEqual(actual, sorted(actual)) + self.assertEqual( + sum((a - b) ** 2 for a, b in zip(values, actual, strict=True)), min(costs) + ) + + def test_invalid(self): + for values in [[True], [1.5]]: + with self.assertRaises(ValueError): + isotonic_regression(values) diff --git a/tests/test_jensen_shannon_divergence.py b/tests/test_jensen_shannon_divergence.py new file mode 100644 index 0000000..bbd50a7 --- /dev/null +++ b/tests/test_jensen_shannon_divergence.py @@ -0,0 +1,27 @@ +import math +import unittest + +from algorithm_lab.jensen_shannon_divergence import jensen_shannon_divergence + + +class Tests(unittest.TestCase): + def test_known_values_and_symmetry(self): + self.assertEqual(jensen_shannon_divergence([1, 0], [0, 1]), 1) + self.assertEqual(jensen_shannon_divergence([1, 3], [2, 6]), 0) + expected = 0.31127812445913283 + self.assertAlmostEqual(jensen_shannon_divergence([1, 0], [1, 1]), expected) + for a, b in [([1, 2, 0], [2, 1, 4]), ([1e308, 0], [1e308, 1e308])]: + self.assertAlmostEqual(jensen_shannon_divergence(a, b), jensen_shannon_divergence(b, a)) + self.assertTrue(math.isfinite(jensen_shannon_divergence([5e-324, 1], [0, 1]))) + + def test_invalid(self): + for a, b in [ + ([], []), + ([1], [1, 2]), + ([0], [1]), + ([-1], [1]), + ([math.nan], [1]), + ([True], [1]), + ]: + with self.assertRaises(ValueError): + jensen_shannon_divergence(a, b) diff --git a/tests/test_kendall_tau.py b/tests/test_kendall_tau.py new file mode 100644 index 0000000..57d0db0 --- /dev/null +++ b/tests/test_kendall_tau.py @@ -0,0 +1,27 @@ +import math +import unittest +from itertools import permutations + +from algorithm_lab.kendall_tau import kendall_tau + + +class Tests(unittest.TestCase): + def test_permutation_inversion_oracle(self): + for n in range(2, 7): + pairs = n * (n - 1) // 2 + for p in permutations(range(n)): + inversions = sum(p[i] > p[j] for i in range(n) for j in range(i + 1, n)) + self.assertAlmostEqual(kendall_tau(range(n), p), (pairs - 2 * inversions) / pairs) + self.assertAlmostEqual(kendall_tau([1, 1, 2], [1, 2, 3]), 2 / math.sqrt(6)) + + def test_invalid(self): + for a, b in [ + ([], []), + ([1], [1]), + ([1, 1], [1, 2]), + ([1, 2], [1]), + ([True, 2], [1, 2]), + ([math.nan, 2], [1, 2]), + ]: + with self.assertRaises(ValueError): + kendall_tau(a, b) diff --git a/tests/test_online_covariance.py b/tests/test_online_covariance.py new file mode 100644 index 0000000..52c0673 --- /dev/null +++ b/tests/test_online_covariance.py @@ -0,0 +1,31 @@ +import math +import random +import statistics +import unittest + +from algorithm_lab.online_covariance import online_covariance + + +class Tests(unittest.TestCase): + def test_standard_library_oracle(self): + rng = random.Random(256) + for n in range(2, 60): + x = [rng.uniform(-100, 100) for _ in range(n)] + y = [rng.uniform(-100, 100) for _ in range(n)] + expected = statistics.covariance(x, y) + self.assertAlmostEqual( + online_covariance(zip(x, y, strict=True), sample=True), expected, places=8 + ) + self.assertAlmostEqual( + online_covariance(zip(x, y, strict=True)), expected * (n - 1) / n, places=8 + ) + self.assertEqual(online_covariance([(3, 5)]), 0) + + def test_invalid(self): + for pairs in [[], [(True, 1)], [(math.nan, 1)], [(1e308, 1e308), (-1e308, -1e308)]]: + with self.assertRaises(ValueError): + online_covariance(pairs) + with self.assertRaises(ValueError): + online_covariance([(1, 2)], sample=True) + with self.assertRaises(ValueError): + online_covariance([(1, 2)], sample=1) diff --git a/tests/test_percentile.py b/tests/test_percentile.py new file mode 100644 index 0000000..6b3c873 --- /dev/null +++ b/tests/test_percentile.py @@ -0,0 +1,30 @@ +import math +import statistics +import unittest + +from algorithm_lab.percentile import percentile + + +class Tests(unittest.TestCase): + def test_standard_library_quantiles(self): + values = [9, -2, 7, 4, 4, 11, 1] + expected = statistics.quantiles(values, n=100, method="inclusive") + for i, value in enumerate(expected, 1): + self.assertAlmostEqual(percentile(values, i / 100), value) + self.assertEqual(percentile([3], 0.4), 3) + self.assertEqual(percentile([-1e308, 1e308], 0.5), 0) + self.assertEqual(percentile(values, 0), -2) + self.assertEqual(percentile(values, 1), 11) + + def test_invalid(self): + for values, q in [ + ([], 0.5), + ([1], -0.1), + ([1], 1.1), + ([math.nan], 0.5), + ([True], 0.5), + ([1], True), + ([10**1000], 0.5), + ]: + with self.assertRaises(ValueError): + percentile(values, q) diff --git a/tests/test_softmax.py b/tests/test_softmax.py new file mode 100644 index 0000000..9a203b6 --- /dev/null +++ b/tests/test_softmax.py @@ -0,0 +1,20 @@ +import math +import unittest + +from algorithm_lab.softmax import softmax + + +class Tests(unittest.TestCase): + def test_probability_and_shift_invariance(self): + for logits in [[-2, 0, 1], [0], [3, 3, 3], [-50, 50]]: + expected = [math.exp(v) / sum(math.exp(x) for x in logits) for v in logits] + for actual, want in zip(softmax(logits), expected, strict=True): + self.assertAlmostEqual(actual, want) + for a, b in zip(softmax(logits), softmax([x + 1000 for x in logits]), strict=True): + self.assertAlmostEqual(a, b) + self.assertEqual(softmax([-1e308, 1e308]), [0, 1]) + + def test_invalid(self): + for values in [[], [math.nan], [math.inf], [True], [10**1000]]: + with self.assertRaises(ValueError): + softmax(values) diff --git a/tests/test_weighted_median.py b/tests/test_weighted_median.py new file mode 100644 index 0000000..0007065 --- /dev/null +++ b/tests/test_weighted_median.py @@ -0,0 +1,29 @@ +import random +import unittest + +from algorithm_lab.weighted_median import weighted_median + + +class Tests(unittest.TestCase): + def test_expanded_sample_oracle(self): + rng = random.Random(251) + for _ in range(100): + values = [rng.randrange(-5, 6) for _ in range(10)] + weights = [rng.randrange(4) for _ in values] + expanded = sorted(v for v, w in zip(values, weights, strict=True) for _ in range(w)) + if expanded: + actual = weighted_median(values, weights) + self.assertEqual(actual, expanded[(len(expanded) - 1) // 2]) + cost = sum(w * abs(v - actual) for v, w in zip(values, weights, strict=True)) + self.assertEqual( + cost, + min( + sum(w * abs(v - c) for v, w in zip(values, weights, strict=True)) + for c in values + ), + ) + + def test_invalid(self): + for a, w in [([], []), ([1], []), ([1], [0]), ([1], [-1]), ([True], [1]), ([1], [True])]: + with self.assertRaises(ValueError): + weighted_median(a, w)