#!/usr/bin/env python3
"""Fit the paper's stated finite-size form to each seeded run.

For each seed folder, read outputs/paper_fss_curves.json (P_inf against p for
N = 150, 500, 2000, 8000, written by the authors' paper_figs_part2.py), find
p_c(N) where P_inf first crosses 0.5 (linear interpolation, as in the authors'
script), then fit p_c(N) = p_c(inf) - c * N**(-3/8) by ordinary least squares
over N = 500, 2000, 8000.

Usage: python3 fit_pc_inf.py            (run from this folder; standard library only)
"""
import json
import pathlib

HERE = pathlib.Path(__file__).resolve().parent
P_GRID = [round(0.15 + 0.02 * i, 3) for i in range(20)]
FIT_SIZES = (500, 2000, 8000)
EXPONENT = 0.375


def p_c_by_size(curves):
    out = {}
    for n in (150, 500, 2000, 8000):
        c = curves[str(n)]
        ys = [c[str(p)] if str(p) in c else c[repr(p)] for p in P_GRID]
        above = [i for i, y in enumerate(ys) if y >= 0.5]
        if above and above[0] > 0:
            i = above[0]
            out[n] = P_GRID[i - 1] + (0.5 - ys[i - 1]) / (ys[i] - ys[i - 1]) * (P_GRID[i] - P_GRID[i - 1])
    return out


def fit(pcs):
    xs = [n ** -EXPONENT for n in FIT_SIZES]
    ys = [pcs[n] for n in FIT_SIZES]
    mx, my = sum(xs) / len(xs), sum(ys) / len(ys)
    slope = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / sum((x - mx) ** 2 for x in xs)
    return my - slope * mx, -slope  # p_c(inf), c


def main():
    print("seed\tp_c(500)\tp_c(2000)\tp_c(8000)\tp_c(inf)\tc")
    for d in sorted(HERE.glob("seed-*")):
        pcs = p_c_by_size(json.loads((d / "outputs" / "paper_fss_curves.json").read_text(encoding="utf-8")))
        pc_inf, c = fit(pcs)
        print(f"{d.name[5:]}\t{pcs[500]:.4f}\t{pcs[2000]:.4f}\t{pcs[8000]:.4f}\t{pc_inf:.4f}\t{c:.4f}")


if __name__ == "__main__":
    main()