The Principia Machina/TPC-PY/gen beta corpus.py: Difference between revisions

Latest revision as of 02:11, 30 September 2026

"""gen_beta_corpus 0.1.0: builds beta_corpus.json for tpc_beta 0.1.0.
Expectations are computed by tpc_beta.py and pinned; hand_checked cases
have a closed form (I_x(1, 1) = x, I_x(2, 1) = x^2, I_x(1, b) = 1 - (1 -
x)^b, I_1/2(a, a) = 1/2, the (1, 1) interval 0.025..0.975, the (2, 1)
quantile sqrt(p), lgamma(1) = lgamma(2) = 0, lgamma(1/2) = ln sqrt(pi))
checked here against the module to 1e-11 relative (the render carries
twelve significant digits) before pinning. Also checks
canon_float against the tree's tpc_strike._canon_float when the tree is
reachable (WORK/tree/tpc)."""
import json
import math
import os
import sys

here = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, here)
from tpc_beta import (betainc, betainc_inv, beta_interval, beta_mean,  # noqa: E402
                      canon_float, lgamma, run_case)

cases = []


def add(name, fn, args, hand=False, check=None):
    if check is not None:
        got = run_case(fn, args)
        v = float(got.split("|")[0])
        assert abs(v - check) <= 1e-11 * max(1.0, abs(check)), (name, got, check)
    cases.append({"name": name, "fn": fn, "args": args, "hand_checked": hand,
                  "expected": run_case(fn, args)})


# ---- mean ------------------------------------------------------------------
add("M01 the prior (1, 1)", "mean", {"a": 1, "b": 1}, True, 0.5)
add("M02 (2, 1)", "mean", {"a": 2, "b": 1}, True, 2 / 3)
add("M03 (11, 3): F 10, A 2 on the prior", "mean", {"a": 11, "b": 3}, True, 11 / 14)
add("M04 (1, 11)", "mean", {"a": 1, "b": 11}, True, 1 / 12)
add("M05 (0.5, 0.5)", "mean", {"a": 0.5, "b": 0.5}, True, 0.5)
add("M06 (100, 1)", "mean", {"a": 100, "b": 1}, True, 100 / 101)
add("M07 (1000000, 1)", "mean", {"a": 1000000, "b": 1}, True)
add("M08 (3.5, 2.25) fractional weights", "mean", {"a": 3.5, "b": 2.25}, True, 3.5 / 5.75)
add("M09 (1.0000001, 1) a float a", "mean", {"a": 1.0000001, "b": 1}, False)
add("M10 (7, 7)", "mean", {"a": 7, "b": 7}, True, 0.5)
add("M11 a = 0 is refused", "mean", {"a": 0, "b": 1}, True)
add("M12 b negative is refused", "mean", {"a": 1, "b": -1}, True)
add("M13 a boolean is refused", "mean", {"a": True, "b": 1}, True)
add("M14 a string is refused", "mean", {"a": "1", "b": 1}, True)
add("M15 (1e-3, 1e-3)", "mean", {"a": 0.001, "b": 0.001}, True, 0.5)
add("M16 (2.5, 0.5)", "mean", {"a": 2.5, "b": 0.5}, True, 2.5 / 3)

# ---- interval --------------------------------------------------------------
add("I01 (1, 1) at 0.95: 0.025 and 0.975", "interval", {"a": 1, "b": 1, "level": 0.95}, True, 0.025)
add("I02 (2, 1) at 0.95: sqrt(0.025) and sqrt(0.975)", "interval", {"a": 2, "b": 1, "level": 0.95}, True, math.sqrt(0.025))
add("I03 (1, 2) at 0.95: 1 - sqrt(0.975) and 1 - sqrt(0.025)", "interval", {"a": 1, "b": 2, "level": 0.95}, True, 1 - math.sqrt(0.975))
add("I04 (11, 3) at 0.95", "interval", {"a": 11, "b": 3, "level": 0.95}, False)
add("I05 (1, 11) at 0.95", "interval", {"a": 1, "b": 11, "level": 0.95}, False)
add("I06 (0.5, 0.5) at 0.95", "interval", {"a": 0.5, "b": 0.5, "level": 0.95}, False)
add("I07 (50, 50) at 0.90", "interval", {"a": 50, "b": 50, "level": 0.9}, False)
add("I08 (100, 1) at 0.99", "interval", {"a": 100, "b": 1, "level": 0.99}, False)
add("I09 (1, 1) at 0.5: 0.25 and 0.75", "interval", {"a": 1, "b": 1, "level": 0.5}, True, 0.25)
add("I10 (3, 1) at 0.95: 0.025^(1/3) and 0.975^(1/3)", "interval", {"a": 3, "b": 1, "level": 0.95}, True, 0.025 ** (1 / 3))
add("I11 level 0 is refused", "interval", {"a": 1, "b": 1, "level": 0}, True)
add("I12 level 1 is refused", "interval", {"a": 1, "b": 1, "level": 1}, True)
add("I13 a = 0 is refused", "interval", {"a": 0, "b": 1, "level": 0.95}, True)
add("I14 (1000, 1000) at 0.95", "interval", {"a": 1000, "b": 1000, "level": 0.95}, False)
add("I15 (1.5, 0.25) at 0.95", "interval", {"a": 1.5, "b": 0.25, "level": 0.95}, False)
add("I16 (26, 4) at 0.95: the thin file's neighbour", "interval", {"a": 26, "b": 4, "level": 0.95}, False)
add("I17 a boolean level is refused", "interval", {"a": 1, "b": 1, "level": True}, True)

# ---- betainc ---------------------------------------------------------------
add("B01 I_x(1, 1) = x at 0.3", "betainc", {"a": 1, "b": 1, "x": 0.3}, True, 0.3)
add("B02 I_x(2, 1) = x^2 at 0.3", "betainc", {"a": 2, "b": 1, "x": 0.3}, True, 0.09)
add("B03 I_x(1, 3) = 1 - (1 - x)^3 at 0.2", "betainc", {"a": 1, "b": 3, "x": 0.2}, True, 1 - 0.8 ** 3)
add("B04 I_1/2(a, a) = 1/2 at a = 7.5", "betainc", {"a": 7.5, "b": 7.5, "x": 0.5}, True, 0.5)
add("B05 x = 0", "betainc", {"a": 2, "b": 3, "x": 0}, True, 0.0)
add("B06 x = 1", "betainc", {"a": 2, "b": 3, "x": 1}, True, 1.0)
add("B07 I_x(1, 1) = x at 0.999", "betainc", {"a": 1, "b": 1, "x": 0.999}, True, 0.999)
add("B08 I_x(5, 1) = x^5 at 0.9", "betainc", {"a": 5, "b": 1, "x": 0.9}, True, 0.9 ** 5)
add("B09 (50.5, 30.25) at 0.6", "betainc", {"a": 50.5, "b": 30.25, "x": 0.6}, False)
add("B10 (0.5, 0.5) at 0.5 (arcsine: 1/2)", "betainc", {"a": 0.5, "b": 0.5, "x": 0.5}, True, 0.5)
add("B11 (1000, 1000) at 0.51", "betainc", {"a": 1000, "b": 1000, "x": 0.51}, False)
add("B12 (1, 1000) at 0.001 = 1 - 0.999^1000", "betainc", {"a": 1, "b": 1000, "x": 0.001}, True, 1 - 0.999 ** 1000)
add("B13 x below 0 is refused", "betainc", {"a": 1, "b": 1, "x": -0.1}, True)
add("B14 x above 1 is refused", "betainc", {"a": 1, "b": 1, "x": 1.1}, True)
add("B15 b = 0 is refused", "betainc", {"a": 1, "b": 0, "x": 0.5}, True)
add("B16 (3, 2) at 0.4 = 6 x^3 (1/3 - x/4) ... closed form 0.1792", "betainc", {"a": 3, "b": 2, "x": 0.4}, True, 4 * 0.4 ** 3 - 3 * 0.4 ** 4)
add("B17 (2, 2) at 0.25 = 3x^2 - 2x^3", "betainc", {"a": 2, "b": 2, "x": 0.25}, True, 3 * 0.25 ** 2 - 2 * 0.25 ** 3)

# ---- betainc_inv -----------------------------------------------------------
add("Q01 p = 0", "betainc_inv", {"a": 2, "b": 3, "p": 0}, True, 0.0)
add("Q02 p = 1", "betainc_inv", {"a": 2, "b": 3, "p": 1}, True, 1.0)
add("Q03 (a, a) at 1/2 is 1/2", "betainc_inv", {"a": 4, "b": 4, "p": 0.5}, True, 0.5)
add("Q04 (2, 1): sqrt(p) at 0.36", "betainc_inv", {"a": 2, "b": 1, "p": 0.36}, True, 0.6)
add("Q05 (1, 3): 1 - (1 - p)^(1/3) at 0.488", "betainc_inv", {"a": 1, "b": 3, "p": 0.488}, True, 0.2)
add("Q06 (1, 1): p itself at 0.025", "betainc_inv", {"a": 1, "b": 1, "p": 0.025}, True, 0.025)
add("Q07 (1, 1): p itself at 0.975", "betainc_inv", {"a": 1, "b": 1, "p": 0.975}, True, 0.975)
add("Q08 (11, 3) at 0.025", "betainc_inv", {"a": 11, "b": 3, "p": 0.025}, False)
add("Q09 (11, 3) at 0.975", "betainc_inv", {"a": 11, "b": 3, "p": 0.975}, False)
add("Q10 (0.5, 0.5) at 0.1", "betainc_inv", {"a": 0.5, "b": 0.5, "p": 0.1}, False)
add("Q11 (1000, 1) at 0.5 = 0.5^(1/1000)", "betainc_inv", {"a": 1000, "b": 1, "p": 0.5}, True, 0.5 ** 0.001)
add("Q12 p below 0 is refused", "betainc_inv", {"a": 1, "b": 1, "p": -0.5}, True)
add("Q13 p above 1 is refused", "betainc_inv", {"a": 1, "b": 1, "p": 1.5}, True)
add("Q14 a negative is refused", "betainc_inv", {"a": -1, "b": 1, "p": 0.5}, True)
add("Q15 (5, 1): p^(1/5) at 0.59049", "betainc_inv", {"a": 5, "b": 1, "p": 0.59049}, True, 0.9)
add("Q16 (30, 70) at 0.5", "betainc_inv", {"a": 30, "b": 70, "p": 0.5}, False)

# ---- lgamma ----------------------------------------------------------------
add("L01 lgamma(1) = 0 (to the approximation's error)", "lgamma", {"x": 1}, True)
add("L02 lgamma(2) = 0 (to the approximation's error)", "lgamma", {"x": 2}, True)
add("L03 lgamma(1/2) = ln sqrt(pi)", "lgamma", {"x": 0.5}, True, math.log(math.sqrt(math.pi)))
add("L04 lgamma(10.5)", "lgamma", {"x": 10.5}, True, math.lgamma(10.5))
add("L05 lgamma(100)", "lgamma", {"x": 100}, True, math.lgamma(100))
add("L06 lgamma(3) = ln 2", "lgamma", {"x": 3}, True, math.log(2))
add("L07 lgamma(6) = ln 120", "lgamma", {"x": 6}, True, math.log(120))
add("L08 lgamma(0.001)", "lgamma", {"x": 0.001}, True, math.lgamma(0.001))
add("L09 lgamma(0.25), reflection", "lgamma", {"x": 0.25}, True, math.lgamma(0.25))
add("L10 lgamma(1e6)", "lgamma", {"x": 1000000}, True, math.lgamma(1000000))
add("L11 lgamma(1.5) = ln(sqrt(pi)/2)", "lgamma", {"x": 1.5}, True, math.log(math.sqrt(math.pi) / 2))
add("L12 x = 0 is refused", "lgamma", {"x": 0}, True)
add("L13 x negative is refused", "lgamma", {"x": -2.5}, True)
add("L14 a boolean is refused", "lgamma", {"x": True}, True)
add("L15 lgamma(171.5)", "lgamma", {"x": 171.5}, True, math.lgamma(171.5))

# ---- canon -----------------------------------------------------------------
tree = os.path.join(here, "..", "..", "tree", "tpc")
strike_canon = None
if os.path.exists(os.path.join(tree, "tpc_strike.py")):
    sys.path.insert(0, tree)
    from tpc_strike import _canon_float as strike_canon  # noqa: E402
cvals = [("K01 zero", 0.0), ("K02 negative zero", -0.0), ("K03 one", 1.0), ("K04 0.1", 0.1), ("K05 1e-5 (exponent form)", 1e-5),
         ("K06 1e12 (exponent form)", 1e12), ("K07 123456789012 (fixed)", 123456789012.0), ("K08 1/3", 1 / 3),
         ("K09 2e-9 (C writes 2e-09)", 2e-9), ("K10 1e22", 1e22), ("K11 0.0001 (fixed)", 0.0001), ("K12 -1.5e-10", -1.5e-10),
         ("K13 an integer argument 7", 7), ("K14 0.025 rounded from 0.024999999999999988", 0.024999999999999988),
         ("K15 999999999999.5 rounds up to 1e12", 999999999999.5), ("K16 -8.881784197001252e-16", -8.881784197001252e-16),
         ("K17 12345.678901234567", 12345.678901234567)]
for name, v in cvals:
    if strike_canon is not None and not isinstance(v, bool) and isinstance(v, float):
        assert canon_float(v) == strike_canon(v), (name, canon_float(v), strike_canon(v))
    add(name, "canon", {"x": v}, True)
add("K18 a boolean is refused", "canon", {"x": True}, True)
if strike_canon is not None:
    import random
    random.seed(20260926)
    same = 0
    for _ in range(100000):
        v = random.uniform(-1, 1) * 10 ** random.randint(-30, 30)
        if random.random() < 0.1:
            v = round(v, random.randint(0, 6))
        if canon_float(v) == strike_canon(v):
            same += 1
    print("canon_float against tpc_strike._canon_float: 17 pinned values and %d/100000 random values identical" % same)
else:
    print("canon_float against tpc_strike._canon_float: tree not reachable, not checked")

json.dump({"corpus": "beta", "version": "0.1.0", "cases": cases},
          open(os.path.join(here, "beta_corpus.json"), "w", encoding="utf-8", newline="\n"), indent=1, sort_keys=True)
print("cases", len(cases))