import math, json, sys
from statistics import NormalDist
T = {"principle": (34, 2157, 27, 1806), "activ400": (10, 817, 9, 774), "activ600": (7, 708, 2, 722), "covidout": (4, 374, 5, 356)}
def tq(p, df):  # Student t quantile via bisection on regularized incomplete beta (simple numeric integration)
    f = lambda u: math.exp(math.lgamma((df+1)/2)-math.lgamma(df/2))/math.sqrt(df*math.pi)*(1+u*u/df)**(-(df+1)/2)
    def cdf(x):
        N = 20000; h = x/N; return 0.5 + sum(f((i+.5)*h) for i in range(N))*h
    lo, hi = 0, 100
    for _ in range(80):
        m = (lo+hi)/2
        if cdf(m) < p: lo = m
        else: hi = m
    return m
def pool(trials):
    y = [math.log((a/n)/(c/m)) for a, n, c, m in trials]; v = [1/a-1/n+1/c-1/m for a, n, c, m in trials]
    w = [1/x for x in v]; k = len(y); fe = sum(wi*yi for wi, yi in zip(w, y))/sum(w)
    Q = sum(wi*(yi-fe)**2 for wi, yi in zip(w, y)); t2 = max(0, (Q-k+1)/(sum(w)-sum(x*x for x in w)/sum(w)))
    ws = [1/(x+t2) for x in v]; mu = sum(wi*yi for wi, yi in zip(ws, y))/sum(ws)
    s = max(1, sum(wi*(yi-mu)**2 for wi, yi in zip(ws, y))/(k-1)); half = tq(0.975, k-1)*math.sqrt(s/sum(ws))
    fh = 1.959963984540054/math.sqrt(sum(w))
    return dict(rr=math.exp(mu), lo=math.exp(mu-half), hi=math.exp(mu+half), tau2=t2, fe=math.exp(fe), fe_lo=math.exp(fe-fh), fe_hi=math.exp(fe+fh), wt_cov=w[-1]/sum(w) if len(w)==4 else None)
base = list(T.values())
r = pool(base); n = sum(a[1]+a[3] for a in base)
print("primary: n", n, "events", sum(a[0] for a in base), sum(a[2] for a in base), {k: round(v, 4) for k, v in r.items() if v is not None})
print("previous 3-contrast:", {k: round(v, 4) for k, v in pool(base[:3]).items() if v is not None})
# TOGETHER unions
lo_env, hi_env, rrs, incl, cnt = 9, 0, [], 0, 0
for ut in range(max(79, 21), min(679, 100)+1):
    for uc in range(max(95, 24), min(679, 119)+1):
        q = pool(base + [(ut, 679, uc, 679)]); cnt += 1; rrs.append(q["rr"])
        lo_env = min(lo_env, q["lo"]); hi_env = max(hi_env, q["hi"]); incl += q["lo"] <= 1 <= q["hi"]
print("TOGETHER scenarios", cnt, "include null", incl, "rr range", round(min(rrs), 4), round(max(rrs), 4), "envelope", round(lo_env, 4), round(hi_env, 4))
# COVID-OUT denominator assignments
rr = [((4+x)/410)/((5+z)/398) for x in range(37) for z in range(43)]
print("COVID-OUT assignments", len(rr), "rr range", round(min(rr), 3), round(max(rr), 3), "all nonevents", round((4/410)/(5/398), 4),
      "first x reversing (z=0)", next(x for x in range(37) if (4+x)/410 > 5/398))