#!/usr/bin/env python3
"""C274: four finite affine systems; exact arithmetic, no parameter search.

Run from anywhere: python3 check_c274.py --out /absolute/output/directory
The frozen packet root is the parent of this script's evidence directory.
"""
import argparse
import csv
import hashlib
import io
import json
from fractions import Fraction as Q
from pathlib import Path

ROOT = Path(__file__).resolve().parent.parent
checks = []


def check(label, condition):
    if not condition:
        raise AssertionError(label)
    checks.append(label)


def read(name):
    return json.loads((ROOT / name).read_text())


def sha(data):
    return hashlib.sha256(data).hexdigest()


def row_reduce(rows):
    a = [[Q(x) for x in r] for r in rows]
    pivots = []
    for col in range(3):
        start = len(pivots)
        pivot = next((r for r in range(start, len(a)) if a[r][col]), None)
        if pivot is None:
            continue
        a[start], a[pivot] = a[pivot], a[start]
        scale = a[start][col]
        a[start] = [x / scale for x in a[start]]
        for r in range(len(a)):
            if r != start:
                scale = a[r][col]
                a[r] = [x - scale * y for x, y in zip(a[r], a[start])]
        pivots.append(col)
    return a, pivots


# A scalar affine expression is (constant, free-parameter coefficient).
def aff(b=0, k=0):
    return Q(b), Q(k)


def add(a, b):
    return a[0] + b[0], a[1] + b[1]


def scale(n, a):
    return Q(n) * a[0], Q(n) * a[1]


def sub(a, b):
    return add(a, scale(-1, b))


def expr(a, parameter):
    b, k = a
    if not k:
        return str(b)
    term = parameter if abs(k) == 1 else f"{abs(k)}*{parameter}"
    if not b:
        return term if k > 0 else "-" + term
    return f"{b}{'+' if k > 0 else '-'}{term}"


def record(a, parameter):
    return {"constant": str(a[0]), "parameter_coefficient": str(a[1]),
            "expression": expr(a, parameter), "determined": a[1] == 0}


def radius(a, parameter):
    return str(abs(a[0])) if a[1] == 0 else f"|{expr(a, parameter)}|"


def main(out):
    scope = read("SCOPE.json")
    manifest = read("SOURCE_MANIFEST.json")
    for name, meta in manifest.items():
        data = (ROOT / "sources" / name).read_bytes()
        check("input integrity: " + name,
              len(data) == meta["bytes"] and sha(data) == meta["sha256"])
    cert = read("sources/C273_CERTIFICATE.json")
    check("C273 predecessor identity",
          manifest["C273_CERTIFICATE.json"]["sha256"] ==
          "c6217745a9340e2380f0a4ec53e9f8e4ebe3b2b16090de2120b00927d92abc23")
    for key, filename in [
        ("report_sha256", "C273_REPORT.md"),
        ("checkpoint_sha256", "C273_CHECKPOINT.md"),
        ("results_sha256", "C273_RESULTS.json"),
        ("source_manifest_sha256", "C273_SOURCE_MANIFEST.json"),
    ]:
        check("C273 certificate binds " + filename,
              cert[key] == manifest[filename]["sha256"])
    prior = read("sources/C273_SOURCE_MANIFEST.json")
    for name in manifest:
        if name in prior:
            check("retained C273 input digest: " + name,
                  prior[name]["sha256"] == manifest[name]["sha256"])

    controls = scope["fixed_controls"]
    source_gap = -Q(controls["Shem_source_offset"])
    increment = Q(controls["Shem_translation_increment"])
    residual = source_gap - increment
    T = Q(controls["retained_source_translation_T"])
    check("fixed 600/598 residual and correspondence",
          (source_gap, increment, residual, T) == (600, 598, 2, 8993))
    endpoints = read("sources/C266_ENDPOINTS.json")
    fields = controls["source_Flood_center_bounds"]
    for profile, bounds in fields.items():
        f = sorted(-r["source_BC_base"] for r in endpoints
                   if r["profile"] == profile and r["source_role"] == "Flood")
        h = sorted(-r["source_BC_base"] for r in endpoints
                   if r["profile"] == profile and r["source_role"] == "Shem"
                   and not r["reflection_after_translation"])
        check(profile + " source field retained",
              f == list(range(bounds[0], bounds[1] + 1)))
        check(profile + " source Shem field 600 earlier",
              h == [x - 600 for x in f])

    m0, m1 = map(Q, fields["MT_ON"])
    l0, l1 = map(Q, fields["LXX_OFF"])
    s0, s1 = map(Q, fields["SP_OFF"])
    equations = scope["premises"]
    # A follows from translated Flood close = reflected translated Shem close.
    check("derive A from source closing coordinates",
          equations[0]["row"] == [2, 0, 0] and
          equations[0]["rhs"] == residual - 2 * m1)
    check("derive B from common shift",
          equations[1]["row"] == [-1, 1, 0] and equations[1]["rhs"] == 0)
    check("derive C from both globally reflected interval endpoints",
          equations[2]["row"] == [1, 0, 1] and
          equations[2]["rhs"] == -(m1 + s0) == -(m0 + s1))
    matrix = [e["row"] + [e["rhs"]] for e in equations]
    reduced, pivots = row_reduce(matrix)
    check("baseline rank three", pivots == [0, 1, 2])
    baseline = [reduced[i][3] for i in range(3)]
    check("unique inherited baseline", baseline == [5290, 5290, 4715])
    determinant = (
        matrix[0][0] * (matrix[1][1]*matrix[2][2] - matrix[1][2]*matrix[2][1])
        - matrix[0][1] * (matrix[1][0]*matrix[2][2] - matrix[1][2]*matrix[2][0])
        + matrix[0][2] * (matrix[1][0]*matrix[2][1] - matrix[1][1]*matrix[2][0])
    )
    check("placement determinant two", determinant == 2)

    families = read("sources/C271_FAMILIES.json")
    coverage = read("sources/C271_COVERAGE.json")
    prior_results = read("sources/C271_RESULTS.json")
    check("retained six source bindings and thirty exclusions",
          (prior_results["cross_pairs"], prior_results["exact_J1_pairs"],
           prior_results["nonmatching_pairs"]) == (36, 6, 30))
    check("retained cumulative coverage",
          all(r["covered_indices"] == list(range(2, 8)) and
              [v["cumulative_source"]["source_index"] for v in r["missing"]] == [0, 1]
              for r in coverage))
    check("six inherited family indices",
          [f["cumulative_index"] for f in families] == list(range(2, 8)))
    for f in families:
        j = f["cumulative_index"]
        check("fixed source family j=" + str(j),
              f["cumulative_Flood_BC_base"] == 5290 - j and
              f["cumulative_Shem_BC_base"] == 5890 - j and
              f["cumulative_Shem_BC_base"] - f["cumulative_Flood_BC_base"] == 600)

    solutions = []
    field_ledger = []
    radius_ledger = []
    regular_ledger = []
    matrix_rows = []
    for case in scope["cases"]:
        cid, omit, parameter = case["id"], case["omit"], case["parameter"]
        selected = [e for e in equations if e["id"] != omit]
        rref, pivots = row_reduce([e["row"] + [e["rhs"]] for e in selected])
        free_columns = [i for i in range(3) if i not in pivots]
        check(cid + " rank and dimension",
              len(pivots) == (3 if omit is None else 2) and
              len(free_columns) == (0 if omit is None else 1))
        direction = [Q(0)] * 3
        if free_columns:
            free = free_columns[0]
            direction[free] = Q(1)
            for r, col in enumerate(pivots):
                direction[col] = -rref[r][free]
            # Orient the first nonzero coefficient to +1, for readable families.
            factor = next(v for v in direction if v)
            direction = [v / factor for v in direction]
        translations = [aff(b, d) for b, d in zip(baseline, direction)]
        residuals = {}
        for equation in equations:
            val = aff(-equation["rhs"])
            for c, trans in zip(equation["row"], translations):
                val = add(val, scale(c, trans))
            residuals[equation["id"]] = record(val, parameter)
            if equation["id"] != omit:
                check(cid + " retained equation " + equation["id"], val == aff())
            else:
                check(cid + " omitted equation not implied", val[0] == 0 and val[1] != 0)
        check(cid + " integral directions preserve annual centers",
              all(x.denominator == 1 for x in baseline + direction) and
              (not free_columns or Q(1) in direction))
        solutions.append({
            "case": cid, "omitted": omit, "rank": len(pivots),
            "degrees_of_freedom": len(free_columns), "parameter": parameter,
            "parameter_domains": {"formal_affine": "Q", "annual_center_admissible": "Z"},
            "rref_augmented": [[str(x) for x in row] for row in rref],
            "reference_origin": [str(x) for x in baseline],
            "null_direction": [str(x) for x in direction],
            "Flood_translations": {k: record(v, parameter) for k, v in zip(["m", "l", "s"], translations)},
            "Shem_translations": {k: record(add(v, aff(increment)), parameter)
                                  for k, v in zip(["m", "l", "s"], translations)},
            "equation_residuals": residuals,
            "free_parameter_selected": False,
        })
        local = {}
        for profile, trans in zip(["MT_ON", "LXX_OFF", "SP_OFF"], translations):
            low, high = map(Q, fields[profile])
            f = [add(trans, aff(low)), add(trans, aff(high))]
            h = [add(trans, aff(low - source_gap + increment)),
                 add(trans, aff(high - source_gap + increment))]
            rh = [scale(-1, h[1]), scale(-1, h[0])]
            local[profile] = f
            check(cid + " " + profile + " field widths and J1",
                  sub(f[1], f[0]) == aff(7) and sub(rh[1], rh[0]) == aff(7) and
                  h == [sub(x, aff(2)) for x in f] and
                  rh == [sub(aff(2), f[1]), sub(aff(2), f[0])])
            # Reflection reverses ancestry/order of the chronological phase pair.
            x = f[1]
            check(cid + " " + profile + " physical phase reflection",
                  [sub(aff(2), add(x, aff(Q(1, 4)))),
                   sub(aff(2), sub(x, aff(Q(1, 4))))] ==
                  [sub(sub(aff(2), x), aff(Q(1, 4))),
                   add(sub(aff(2), x), aff(Q(1, 4)))])
            field_ledger.append({
                "case": cid, "profile": profile, "parameter": parameter,
                "Flood_centers": [record(x, parameter) for x in f],
                "unreflected_Shem_centers": [record(x, parameter) for x in h],
                "reflected_Shem_centers_in_chronological_order": [record(x, parameter) for x in rh],
                "field_center_width": "7", "individual_seasonal_envelope_length": "15/2",
                "local_reflection": "J1(x)=2-x",
                "full_field_radius_about_1": "max(" + ",".join(radius(sub(x, aff(1)), parameter) for x in f) + ")",
            })
        m, l, s = translations
        mt, lxx, sp = (local[x] for x in ["MT_ON", "LXX_OFF", "SP_OFF"])
        join_gap = sub(sub(aff(2), mt[1]), mt[1])
        cross_offset = sub(lxx[0], mt[0])
        mirror_error = add(sp[0], mt[1])
        check(cid + " MT join gap from source equation",
              join_gap == sub(aff(10580), scale(2, m)))
        check(cid + " common-shift corresponding endpoint difference",
              cross_offset == add(aff(6), sub(l, m)) and
              sub(lxx[1], mt[1]) == cross_offset)
        check(cid + " both SP global-reflection residuals agree",
              mirror_error == add(sp[1], mt[0]) ==
              sub(add(m, s), aff(10005)))
        selected_radii = {}
        for label, x in [("MT_opening", mt[0]), ("SP_closing", sp[1]), ("LXX_closing", lxx[1])]:
            displacement = sub(x, aff(1))
            selected_radii[label] = {
                "source_selected_center": record(x, parameter),
                "signed_displacement_from_1": record(displacement, parameter),
                "radius": radius(displacement, parameter),
                "radius_determined": displacement[1] == 0,
                "selected_endpoint_guaranteed_outermost_for_all_parameters": displacement[1] == 0,
            }
        radius_ledger.append({
            "case": cid, "parameter": parameter,
            "inherited_selected_endpoints": selected_radii,
            "MT_closing_join_signed_gap": record(join_gap, parameter),
            "LXX_minus_MT_corresponding_Flood_centers": record(cross_offset, parameter),
            "SP_global_Mirror_of_MT_Flood_endpoint_error": record(mirror_error, parameter),
        })

        # Compose retained T^-1 with the cumulative LXX maps.
        noah_intercept = sub(l, aff(T))
        flood_intercept = sub(aff(T - increment), l)
        delta = sub(l, aff(5290))
        check(cid + " universal regular binding residual",
              sub(add(noah_intercept, flood_intercept), aff(2)) == aff(-600))
        case_families = []
        for f in families:
            j = f["cumulative_index"]
            noah_source = T - f["cumulative_Flood_BC_base"]
            flood_source = T - f["cumulative_Shem_BC_base"]
            n = add(aff(noah_source), noah_intercept)
            z = add(aff(-flood_source), flood_intercept)
            check(cid + " retained source pair j=" + str(j),
                  noah_source - flood_source == 600 and
                  n == add(aff(j), delta) and z == sub(aff(2), n))
            case_families.append({
                "cumulative_index": j, "gear": f["gear"], "binding": f["binding"],
                "Noah_source_center": str(noah_source), "Flood_source_center": str(flood_source),
                "local_Noah_center": record(n, parameter),
                "local_Flood_center": record(z, parameter),
                "signed_Noah_displacement_from_1": record(sub(n, aff(1)), parameter),
                "radius": radius(sub(n, aff(1)), parameter),
            })
        if cid == "baseline":
            original_maps = read("sources/C270_MAPS.json")
            check("recover original induced maps",
                  original_maps["regular_reflected_Noah_map"] == ["1", str(noah_intercept[0])] and
                  original_maps["regular_reflected_Flood_map"] == ["-1", str(flood_intercept[0])])
        signed_close = sub(lxx[1], aff(1))
        eta_components = [scale(1000, add(signed_close, aff(d))) for d in [-Q(1, 4), Q(1, 4)]]
        regular_ledger.append({
            "case": cid, "parameter": parameter, "delta_l": record(delta, parameter),
            "source_T": "q+8993",
            "regular_Noah_map": {"slope": "1", "intercept": record(noah_intercept, parameter)},
            "regular_Flood_map": {"slope": "-1", "intercept": record(flood_intercept, parameter)},
            "Noah_subtracted_shift": record(scale(-1, noah_intercept), parameter),
            "all_pair_residual": "RN(yN)+RF(yF)-2=yN-yF-600",
            "source_match_count_retained_without_new_scan": 6,
            "source_nonmatch_count_retained_without_new_scan": 30,
            "cumulative_indices_retained": list(range(2, 8)),
            "missing_source_indices_retained": [0, 1],
            "matched_families": case_families,
            "source_selected_closing_Noah_side_signed_eta_components":
                [record(x, parameter) for x in eta_components],
            "C271_1_to_6_radii_fixed_in_inherited_order": delta[1] == 0,
        })
        matrix_rows.append({
            "case": cid, "rank": len(pivots), "free_parameters": len(free_columns),
            "MT_translation_fixed": m[1] == 0, "LXX_translation_fixed": l[1] == 0,
            "cumulative_SP_translation_fixed": s[1] == 0,
            "local_J1_all_profiles": True,
            "MT_closing_join_for_all_values": join_gap == aff(),
            "SP_global_Mirror_of_MT_Flood_for_all_values": mirror_error == aff(),
            "MT_LXX_corresponding_offset_6_for_all_values": cross_offset == aff(6),
            "MT_selected_radius_7_fixed": selected_radii["MT_opening"]["radius_determined"],
            "SP_selected_radius_5_fixed": selected_radii["SP_closing"]["radius_determined"],
            "LXX_selected_radius_6_fixed": selected_radii["LXX_closing"]["radius_determined"],
            "regular_source_six_bindings_fixed": True,
            "regular_local_radii_1_to_6_fixed": delta[1] == 0,
            "regular_Noah_shift_7K_fixed": noah_intercept == aff(-3703),
        })

    check("case directions recovered",
          [r["null_direction"] for r in solutions] ==
          [["0", "0", "0"], ["1", "1", "-1"], ["0", "1", "0"], ["0", "0", "1"]])
    check("all selected radius statuses",
          [[r[x] for x in ["MT_selected_radius_7_fixed", "SP_selected_radius_5_fixed",
                          "LXX_selected_radius_6_fixed"]] for r in matrix_rows] ==
          [[True, True, True], [False, False, False], [True, True, False], [True, False, True]])
    check("cumulative SP omission leaves regular SP placement fixed",
          matrix_rows[3]["regular_local_radii_1_to_6_fixed"] and
          not matrix_rows[3]["cumulative_SP_translation_fixed"])
    outputs = {
        "EQUATION_SOLUTIONS.json": {"variable_order": ["m", "l", "s"],
                                    "determinant": determinant, "cases": solutions},
        "FIELD_AND_SYMMETRY_LEDGER.json": field_ledger,
        "RADIUS_LEDGER.json": radius_ledger,
        "REGULAR_MAPS.json": regular_ledger,
        "RESULTS.json": {
            "step": "C274", "status": "passed", "check_count": len(checks), "checks": checks,
            "frozen_inputs": len(manifest), "case_count": len(solutions),
            "case_ranks": [r["rank"] for r in solutions],
            "case_degrees_of_freedom": [r["degrees_of_freedom"] for r in solutions],
            "placement_premises_irredundant_within_held_affine_system": True,
            "historical_necessity_of_placement_premises_proved": False,
            "fixed_598_increment_audited_or_removed": False,
            "free_parameters_selected_or_fitted": False,
            "new_source_records": 0, "historical_campaigns_rerun": 0,
            "regular_pair_invariance_method": "universal residual identity, using retained C271 counts; no new 36-pair scan",
            "new_independent_review": False, "next_step_executed": False,
        },
    }
    out.mkdir(parents=True, exist_ok=True)
    for name, obj in outputs.items():
        (out / name).write_text(json.dumps(obj, ensure_ascii=False, indent=2) + "\n")
    buf = io.StringIO(newline="")
    writer = csv.DictWriter(buf, fieldnames=list(matrix_rows[0]))
    writer.writeheader()
    writer.writerows(matrix_rows)
    (out / "INVARIANT_MATRIX.csv").write_text(buf.getvalue())
    print(json.dumps({"status": "passed", "checks": len(checks),
                      "cases": len(solutions), "output_files": len(outputs) + 1}))


if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument("--out", type=Path, required=True)
    main(parser.parse_args().out)