# Source: content/notes/math/numerical-methods.md
# Independent CPU example; use the curriculum environment.
# See /notes/ml/#example-environment or /notes/deep-learning/#example-environment.

import math
import numpy as np
from scipy.special import logsumexp

tiny = np.finfo(np.float16).smallest_subnormal
assert tiny == np.float16(2.**-24)
assert np.float16(1e-5) > 0  # subnormal is representable on this CPU path
assert np.spacing(np.float32(1e8)) == 8.
lost = np.array([1e8, 1e8+1, 1e8+2], dtype=np.float32)
assert np.var(lost.astype(np.float64)) == 0.
kept = np.array([10000., 10001., 10002.], dtype=np.float32)
assert np.isclose(np.var(kept.astype(np.float64)), 2/3)

def pairwise(values):
    if len(values) <= 1:
        return sum(values)
    middle = len(values)//2
    return pairwise(values[:middle])+pairwise(values[middle:])

def neumaier(values):
    total, correction = 0., 0.
    for x in values:
        new = total+x
        if abs(total) >= abs(x):
            correction += (total-new)+x
        else:
            correction += (x-new)+total
        total = new
    return total+correction

values = [1e16, 1., -1e16]
naive = 0.
for value in values:
    naive += value
reference = math.fsum(values)
assert reference == 1. and naive == 0.
assert neumaier(values) == reference
print("naive / pairwise / compensated / fsum:",
      naive, pairwise(values), neumaier(values), reference)
assert np.isneginf(logsumexp([-np.inf, -np.inf]))
assert np.isposinf(logsumexp([np.inf, 0.]))
# Older SciPy versions reject empty reductions; apply the explicit sum convention.
empty = np.array([], dtype=float)
empty_lse = -np.inf if empty.size == 0 else logsumexp(empty)
assert np.isneginf(empty_lse)