"""Independent recomputation of the age-35/40 induction RD from data/cells_main.csv (pooled over years)."""
import csv, sys, math, collections
rows = list(csv.DictReader(open(sys.argv[1] + "/data/cells_main.csv")))
agg = collections.defaultdict(lambda: [0.0, 0.0])
for r in rows:
    if r["nullip"] == "1":
        a = int(r["age"]); agg[a][0] += float(r["ind_n"]); agg[a][1] += float(r["ind_k"])
rate = {a: v[1] / v[0] for a, v in agg.items() if v[0] > 0}
print("nulliparous induction % by age:", {a: round(100*rate[a], 2) for a in range(28, 45)})
print("second differences (pp/yr^2):", {a: round(100*(rate[a+1]-2*rate[a]+rate[a-1]), 2) for a in range(27, 44)})
def rd(c, h=3, M=0.001):
    allw = []
    for sign, rng in ((-1, range(c-h, c)), (1, range(c, c+h))):
        pts = [(a + 0.5 - c, rate[a], agg[a][0]) for a in rng]
        ker = [(x, y, n, n * max(0, 1 - abs(x)/h)) for x, y, n in pts]
        sw = sum(k[3] for k in ker); sx = sum(k[3]*k[0] for k in ker)/sw
        sxx = sum(k[3]*(k[0]-sx)**2 for k in ker)
        for x, y, n, w in ker:
            allw.append((x, sign * (w/sw - w*(x - sx)*sx/sxx), y, n))
    tau = sum(l*y for x, l, y, n in allw)
    se = math.sqrt(sum(l*l*y*(1-y)/n for x, l, y, n in allw))
    bias = M/2*sum(abs(l)*x*x for x, l, y, n in allw)
    return tau, se, bias
for c in (35, 40):
    for M in (0.001, 0.002, 0.004):
        t, s, b = rd(c, M=M)
        print(f"c={c} M={M*100:.1f}pp: tau {100*t:.2f} se {100*s:.3f} bias {100*b:.2f} CI {100*(t-b-1.96*s):.2f}..{100*(t+b+1.96*s):.2f}")
    t, s, b = rd(c); print(f"  M at which the age-{c} honest CI reaches 0: {(t-1.96*s)/(b/0.1):.2f} pp/yr^2")
print("placebo cutoffs at M=0.1pp:")
for c in (31, 32, 33, 37, 38):
    t, s, b = rd(c); print(f"  c={c}: tau {100*t:.2f} CI {100*(t-b-1.96*s):.2f}..{100*(t+b+1.96*s):.2f}")
for a in (33, 38, 42):
    se2 = math.sqrt(sum(rate[x]*(1-rate[x])/agg[x][0]*m*m for x, m in ((a-1,1),(a,2),(a+1,1))))
    print(f"second difference at {a}: {100*(rate[a+1]-2*rate[a]+rate[a-1]):.2f} pp, sampling SE {100*se2:.2f} pp")