#!/usr/bin/env python3
"""Tier-1 item 2 -- the homogeneous-core control.

Why this run matters
--------------------
The reported endpoint is an area-weighted mean over the OUTER surface. Core
heterogeneity plausibly averages out on its way there, and Main Fig. 3c already
hints at exactly that: the Spearman rho between the chromatin modulus scale and
the paired change is only -0.08, i.e. scaling every chromatin modulus by a
factor of four barely moves the global endpoint. A reviewer will ask whether the
digitised Reynolds map does any causal work for the headline claim.

Both outcomes help the paper:
  * endpoint unchanged -> say so plainly; heterogeneity drives the local fields
    (already exploratory) and the paper is a natural-volume-shrinkage study;
  * endpoint moves -> that is the actual finding and belongs in Main Fig. 2.

Implementation
--------------
The solver already has a `homogeneous` branch (`core_heterogeneity_model`),
but the runner hard-codes `REYNOLDS_2021_MODEL`, so this script patches
`_options` rather than editing upstream source.

Choosing the equivalent modulus is the part that decides whether the control is
fair. `--core-young-pa 800` (the upstream fallback, and the source of the
Core-Young/Poisson entries in Table S1) is NOT a fair equivalent: it is an
unrelated legacy value. Run the Voigt and Reuss equivalents instead -- they
bracket any physically admissible homogenisation of the same group population,
so if the endpoint sits inside a narrow bracket the conclusion is safe either
way. Use 04_equivalent_core_modulus.py to compute them, or pass them directly.
"""

from __future__ import annotations

import argparse
import json
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parent))
from _common import (  # noqa: E402
    add_common_arguments,
    import_runner,
    paired_change,
    patch_options,
    read_endpoints,
)


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__)
    add_common_arguments(parser)
    parser.add_argument(
        "--core-young-pa",
        type=float,
        nargs="+",
        required=True,
        metavar="E",
        help=(
            "One or more homogeneous core Young moduli (Pa). Supply the Voigt "
            "and Reuss equivalents from 04_equivalent_core_modulus.py. Adding "
            "800 reproduces the upstream legacy value for comparison only."
        ),
    )
    parser.add_argument(
        "--label",
        nargs="*",
        default=None,
        help="Optional labels matching --core-young-pa (default: E<value>).",
    )
    parser.add_argument(
        "--element-order",
        choices=("p1", "p2"),
        default="p2",
    )
    args = parser.parse_args()

    labels = args.label or [f"E{value:g}" for value in args.core_young_pa]
    if len(labels) != len(args.core_young_pa):
        raise SystemExit("--label must have the same count as --core-young-pa")

    print("Homogeneous-core control -- planned runs")
    print(f"  project root : {args.project_root}")
    print(f"  output root  : {args.out_root}")
    print(f"  primary mesh : {args.primary_mesh}")
    for label, modulus in zip(labels, args.core_young_pa, strict=True):
        print(f"  homogeneous_{label:<12} core E = {modulus:g} Pa")
    print(
        "  Poisson ratio comes from the upstream core config; keep it equal to "
        "the heterogeneous runs' effective value (K/G = 30 -> nu = 0.4835), not "
        "the 0.45 printed in Table S1."
    )
    if args.dry_run:
        print("\n--dry-run: nothing solved.")
        return 0

    runner = import_runner(args.project_root)
    patch_options(runner, core_heterogeneity_model="homogeneous")
    args.out_root.mkdir(parents=True, exist_ok=True)

    discretization = f"{args.primary_mesh}_{'P2' if args.element_order == 'p2' else 'P1'}"
    results = {}
    for label, modulus in zip(labels, args.core_young_pa, strict=True):
        out_dir = args.out_root / f"homogeneous_{label}"
        print(f"\n=== homogeneous core E = {modulus:g} Pa -> {out_dir} ===", flush=True)
        runner.run_analysis(
            out_dir,
            primary_mesh=args.primary_mesh,
            core_young_pa=modulus,
            run_intermediate_p2=(args.element_order == "p2"),
            resume=args.resume,
            sparse_backend=args.sparse_backend,
            rigid_constraint_linear_solver=args.rigid_constraint_linear_solver,
        )
        try:
            rows = read_endpoints(out_dir)
            control, hyper, difference, percent = paired_change(rows, discretization)
            results[label] = {
                "core_young_pa": modulus,
                "control_mn_per_m": control,
                "hyperosmotic_mn_per_m": hyper,
                "absolute_difference_mn_per_m": difference,
                "percent_change": percent,
                "out_dir": str(out_dir),
            }
        except (FileNotFoundError, KeyError) as exc:
            print(f"  warning: could not read {label}: {exc}")

    output = args.out_root / "homogeneous_control.json"
    output.write_text(
        json.dumps(
            {
                "discretization": discretization,
                "results": results,
                "interpretation": (
                    "Compare percent_change against the heterogeneous run at the "
                    "same discretisation. A difference smaller than the "
                    "discretisation spread (>=3.6 % same-mesh P1/P2) means the "
                    "global endpoint cannot distinguish a heterogeneous core "
                    "from its homogeneous equivalent."
                ),
            },
            indent=2,
        ),
        encoding="utf-8",
    )
    print(f"\nWrote {output}")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())