#!/usr/bin/env python
"""Three-region M82 tutorial with an 8-epoch joint center fit.

North and south preserve the accepted regional baselines. The center is not fit
as one merged static spectrum: eight Chandra epochs share diffuse plasma/CX/Fe-K
parameters while each epoch has its own absorbed XRB/ULX continuum. A weighted
mean source is propagated through the SQUDE response.
"""
from __future__ import annotations

import argparse
import copy
import hashlib
import json
import os
import sys
from pathlib import Path
from types import SimpleNamespace
from typing import Any

os.environ.setdefault("MPLCONFIGDIR", "/tmp/m82-three-region-joint-mpl")
os.environ.setdefault("XDG_CACHE_HOME", "/tmp/m82-three-region-joint-cache")

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import scienceplots  # noqa: F401

from sherpa.astro.fake import poisson_noise
from sherpa.astro.ui import (
    calc_energy_flux,
    clean,
    create_model_component,
    fake_pha,
    fit,
    freeze,
    get_data,
    get_data_plot,
    get_fit_results,
    get_model_component,
    get_model_plot,
    group_counts,
    load_pha,
    notice_id,
    save_pha,
    set_analysis,
    set_method,
    set_rng,
    set_source,
    set_stat,
    ungroup,
)

ROOT = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(ROOT / "scripts"))
import make_m82_mock_tutorial_figure_set as core  # noqa: E402
import make_m82_three_region_mock_tutorial_figures as precursor  # noqa: E402

OUT_DIR = ROOT / "data" / "fits" / "m82_three_region_mock_tutorial"
CENTER_ROOT = precursor.SPEC_ROOT / "v40_center"
CENTER_JOINT_PHOTON_INDEX_INIT = 1.6
CENTER_JOINT_NH_INIT = [7.0080, 6.8376, 3.9863, 3.6950, 4.0596, 5.8073, 21.6074, 2.7040]
CENTER_JOINT_AMPL_INIT = [0.0023301, 0.0022742, 0.0043836, 0.0052007, 0.0053063, 0.0033966, 0.0025452, 0.0049955]
CENTER_STYLE = {
    **core.COMPONENT_STYLE,
    "fe64": ("Fe Kα 6.40 keV", "#CC79A7", ":"),
    "xrb": ("epoch-variable XRB/ULX average", "#7F3C8D", "-."),
}


def sha256(path: Path) -> str:
    digest = hashlib.sha256()
    with path.open("rb") as handle:
        for chunk in iter(lambda: handle.read(1024 * 1024), b""):
            digest.update(chunk)
    return digest.hexdigest()


def provenance_path(path: Path, output_dir: Path) -> str:
    try:
        return str(path.relative_to(ROOT))
    except ValueError:
        try:
            return str(path.relative_to(output_dir))
        except ValueError:
            return path.name


def add_terms(terms: list[Any]) -> Any:
    total = terms[0]
    for term in terms[1:]:
        total = total + term
    return total


def setup_center_joint_model() -> tuple[dict[int, Any], dict[str, Any], Any, dict[int, Any]]:
    bundle = core.setup_literature_model(1)
    gal = get_model_component("gal")
    absth = get_model_component("absth")
    vcool = get_model_component("vcool")
    vhot = get_model_component("vhot")
    gce_o = get_model_component("gce_o")
    gce_mg = get_model_component("gce_mg")

    fe64 = create_model_component("xsgaussian", "fe64")
    fe64.LineE = 6.40
    freeze(fe64.LineE)
    fe64.Sigma = 0.001
    freeze(fe64.Sigma)
    fe64.norm = 1.0e-6

    shared = {
        "cool_vapec": gal * absth * vcool,
        "hot_vapec": gal * absth * vhot,
        "cx_o": gal * absth * gce_o,
        "cx_mg": gal * absth * gce_mg,
        "fe64": gal * absth * fe64,
    }
    shared_total = add_terms(list(shared.values()))
    sources: dict[int, Any] = {}
    xrb_terms: dict[int, Any] = {}
    for dataset_id in range(1, 9):
        absorber = create_model_component("xstbabs", f"absxrb{dataset_id}")
        xrb = create_model_component("powlaw1d", f"xrb{dataset_id}")
        absorber.nH = CENTER_JOINT_NH_INIT[dataset_id - 1]
        xrb.ref = 1.0
        freeze(xrb.ref)
        xrb.ampl = CENTER_JOINT_AMPL_INIT[dataset_id - 1]
        if dataset_id == 1:
            xrb.gamma = CENTER_JOINT_PHOTON_INDEX_INIT
            xrb.gamma.min = 0.5
            xrb.gamma.max = 3.0
        else:
            xrb.gamma = get_model_component("xrb1").gamma
        xrb_terms[dataset_id] = gal * absorber * xrb
        sources[dataset_id] = shared_total + xrb_terms[dataset_id]
        set_source(dataset_id, sources[dataset_id])
    return sources, shared, shared_total, xrb_terms


def fold_components(dataset_id: int, total_source: Any, components: dict[str, Any]) -> tuple[Any, dict[str, Any]]:
    set_source(dataset_id, total_source)
    total_plot = copy.deepcopy(get_model_plot(dataset_id))
    plots: dict[str, Any] = {}
    for name, source in components.items():
        set_source(dataset_id, source)
        plots[name] = copy.deepcopy(get_model_plot(dataset_id))
    set_source(dataset_id, total_source)
    return total_plot, plots


def aggregate_plot(plots: list[Any], weights: np.ndarray, target: Any) -> Any:
    values = np.zeros_like(np.asarray(target.x, dtype=float))
    for plot, weight in zip(plots, weights, strict=True):
        values += weight * core.rebin_density_conserving_rate(
            np.asarray(plot.xlo), np.asarray(plot.xhi), np.asarray(plot.y),
            np.asarray(target.xlo), np.asarray(target.xhi),
        )
    return SimpleNamespace(
        x=np.asarray(target.x), xlo=np.asarray(target.xlo), xhi=np.asarray(target.xhi),
        y=values,
    )


def positive_limits(arrays: list[np.ndarray]) -> tuple[float, float]:
    positive = np.concatenate([np.asarray(values)[np.asarray(values) > 0] for values in arrays])
    return max(1.0e-6, float(np.percentile(positive, 1)) * 0.45), float(np.max(positive)) * 3.0


def make_center_chandra_figure(data: Any, model: Any, components: dict[str, Any], path: Path, width: float) -> None:
    model_y = core.model_on_data_grid(model, data)
    component_y = {name: core.model_on_data_grid(plot, data) for name, plot in components.items()}
    ymin, ymax = positive_limits([np.asarray(data.y), model_y, *component_y.values()])
    fig, ax = plt.subplots(figsize=(width, width * 0.68))
    for name, values in component_y.items():
        label, color, linestyle = CENTER_STYLE[name]
        ax.plot(data.x, values, color=color, ls=linestyle, lw=0.75 if width < 4 else 1.0, label=label)
    ax.plot(data.x, model_y, color="#D55E00", lw=1.4, label="8-epoch joint folded total")
    ax.errorbar(data.x, data.y, yerr=data.yerr, fmt=".", ms=2.2, color="0.10", label="combined Chandra ACIS data")
    ax.set_yscale("log")
    ax.set_xlim(0.3, 7.05)
    ax.set_ylim(ymin, ymax)
    ax.set_xlabel("Energy (keV)")
    ax.set_ylabel(r"Counts s$^{-1}$ keV$^{-1}$")
    ax.set_title('M82 Center starburst 40"×80": joint 8-epoch Chandra model', fontsize=7.9 if width < 4 else 10.5, pad=72 if width < 4 else 6)
    if width < 4:
        ax.legend(loc="lower center", bbox_to_anchor=(0.5, 1.01), fontsize=4.0, ncol=2, frameon=True, columnspacing=0.7, handlelength=2.3)
    else:
        ax.legend(loc="upper right", fontsize=6.2, ncol=2, frameon=True)
    ax.grid(alpha=0.20)
    fig.savefig(path, dpi=320, bbox_inches="tight")
    plt.close(fig)


def fit_parameter_dict(result: Any) -> dict[str, float]:
    return {name: float(value) for name, value in zip(result.parnames, result.parvals, strict=True)}


def process_center(region: precursor.RegionSpec, output_dir: Path) -> dict[str, Any]:
    clean()
    epoch_pha = sorted(CENTER_ROOT.glob("v40_center_early100ks_src*.pi"))
    if len(epoch_pha) != 8:
        raise RuntimeError(f"Expected 8 center epoch PHAs, found {len(epoch_pha)}")

    epochs: list[dict[str, Any]] = []
    for dataset_id, pha in enumerate(epoch_pha, 1):
        load_pha(dataset_id, str(pha))
        set_analysis(dataset_id, "energy")
        notice_id(dataset_id, 0.3, 7.0)
        group_counts(dataset_id, 20)
        data = get_data(dataset_id)
        epochs.append({
            "dataset_id": dataset_id,
            "pha": pha,
            "exposure_s": float(data.exposure),
            "raw_counts": int(np.sum(data.counts)),
            "count_rate": float(np.sum(data.counts) / data.exposure),
        })

    set_stat("chi2gehrels")
    set_method("levmar")
    sources, shared_components, shared_total, xrb_terms = setup_center_joint_model()
    fit(*range(1, 9))
    fit_result = get_fit_results()

    exposures = np.asarray([epoch["exposure_s"] for epoch in epochs], dtype=float)
    weights = exposures / np.sum(exposures)
    average_xrb = add_terms([float(weight) * xrb_terms[i] for i, weight in enumerate(weights, 1)])
    average_source = shared_total + average_xrb
    average_components = {**shared_components, "xrb": average_xrb}

    combined_pha = CENTER_ROOT / "v40_center_early100ks_combined_src.pi"
    load_pha(9, str(combined_pha))
    set_analysis(9, "energy")
    notice_id(9, 0.3, 7.0)
    group_counts(9, 20)
    combined_data = copy.deepcopy(get_data_plot(9))
    combined_exposure_s = float(get_data(9).exposure)
    exposure_sum_s = float(np.sum(exposures))
    exposure_sum_relative_error = abs(combined_exposure_s - exposure_sum_s) / max(
        abs(exposure_sum_s), 1.0e-30
    )
    if exposure_sum_relative_error > 1.0e-10:
        raise RuntimeError(
            "Center epoch exposure sum does not match combined PHA: "
            f"{exposure_sum_relative_error:.3e}"
        )

    totals: list[Any] = []
    epoch_explicit_validation: list[dict[str, float]] = []
    folded_by_name: dict[str, list[Any]] = {name: [] for name in average_components}
    for dataset_id in range(1, 9):
        epoch_components = {**shared_components, "xrb": xrb_terms[dataset_id]}
        total_plot, component_plots = fold_components(dataset_id, sources[dataset_id], epoch_components)
        epoch_pha_path = epochs[dataset_id - 1]["pha"]
        epoch_explicit_validation.append(precursor.validate_explicit_response_fold(
            sources[dataset_id],
            total_plot,
            epoch_pha_path.with_suffix(".arf"),
            epoch_pha_path.with_suffix(".rmf"),
        ))
        totals.append(total_plot)
        for name, plot in component_plots.items():
            folded_by_name[name].append(plot)
    combined_model = aggregate_plot(totals, weights, combined_data)
    combined_components = {
        name: aggregate_plot(plots, weights, combined_data)
        for name, plots in folded_by_name.items()
    }
    combined_closure = core.validate_component_closure(combined_model, combined_components)

    make_center_chandra_figure(combined_data, combined_model, combined_components, output_dir / "v40_center_chandra_data_model_components_1col.png", 3.35)
    make_center_chandra_figure(combined_data, combined_model, combined_components, output_dir / "v40_center_chandra_data_model_components_2col.png", 7.00)

    set_rng(np.random.default_rng(region.seed))
    set_source(10, average_source)
    squde_pha = output_dir / "v40_center_squde_50ks.pha"
    fake_pha(10, arf=str(precursor.SQUDE_ARF), rmf=str(precursor.SQUDE_RMF), exposure=precursor.EXPOSURE_S, grouped=False, method=poisson_noise)
    save_pha(10, str(squde_pha), clobber=True)
    raw_channels = np.asarray(get_data(10).channel, dtype=float).copy()
    raw_counts_vector = np.asarray(get_data(10).counts, dtype=float).copy()
    raw_squde_counts = int(np.sum(raw_counts_vector))
    set_analysis(10, "energy")
    notice_id(10, 0.3, 4.0)
    group_counts(10, 50)
    grouping_min50 = np.asarray(get_data(10).grouping, dtype=np.int64).copy()
    squde_data = copy.deepcopy(get_data_plot(10))
    squde_model, squde_components = fold_components(10, average_source, average_components)
    squde_closure = core.validate_component_closure(squde_model, squde_components)
    explicit_squde_fold = precursor.validate_explicit_response_fold(
        average_source,
        squde_model,
        precursor.SQUDE_ARF,
        precursor.SQUDE_RMF,
    )
    squde_kernel = core.read_rmf_kernel(precursor.SQUDE_RMF)
    squde_target_lo = np.asarray(squde_model.xlo, dtype=float)
    squde_target_hi = np.asarray(squde_model.xhi, dtype=float)
    explicit_average = core.rebin_density_conserving_rate(
        squde_kernel.channel_lo,
        squde_kernel.channel_hi,
        core.fold_source_with_kernel(average_source, squde_kernel, precursor.SQUDE_ARF),
        squde_target_lo,
        squde_target_hi,
    )
    weighted_epoch_squde = np.zeros_like(explicit_average)
    for dataset_id, weight in enumerate(weights, 1):
        epoch_density = core.fold_source_with_kernel(
            sources[dataset_id], squde_kernel, precursor.SQUDE_ARF
        )
        weighted_epoch_squde += weight * core.rebin_density_conserving_rate(
            squde_kernel.channel_lo,
            squde_kernel.channel_hi,
            epoch_density,
            squde_target_lo,
            squde_target_hi,
        )
    linearity_relative = np.abs(weighted_epoch_squde - explicit_average) / np.maximum(
        np.abs(explicit_average), 1.0e-30
    )
    squde_linearity = {
        "weight_sum": float(np.sum(weights)),
        "median_relative_error": float(np.median(linearity_relative)),
        "p95_relative_error": float(np.percentile(linearity_relative, 95)),
        "max_relative_error": float(np.max(linearity_relative)),
    }
    if squde_linearity["p95_relative_error"] > 1.0e-12:
        raise RuntimeError(
            "Center exposure-weighted SQUDE source linearity failed: "
            f"p95={squde_linearity['p95_relative_error']:.3e}"
        )

    epoch_controls: list[np.ndarray] = []
    epoch_control_validation: list[dict[str, Any]] = []
    for dataset_id, epoch in enumerate(epochs, 1):
        epoch_rmf = epoch["pha"].with_suffix(".rmf")
        control, validation = precursor.build_chandra_rmf_control(
            sources[dataset_id],
            epoch_rmf,
            squde_model,
        )
        epoch_controls.append(control)
        epoch_control_validation.append({
            "dataset_id": dataset_id,
            "exposure_weight": float(weights[dataset_id - 1]),
            "rmf": precursor.provenance_path(epoch_rmf),
            "rmf_sha256": sha256(epoch_rmf),
            **validation,
        })
    chandra_control = np.sum(
        np.asarray(weights)[:, None] * np.asarray(epoch_controls),
        axis=0,
    )
    target_width = np.asarray(squde_model.xhi, dtype=float) - np.asarray(squde_model.xlo, dtype=float)
    weighted_native_rate = float(sum(
        weight * validation["native_overlap_rate_counts_s"]
        for weight, validation in zip(weights, epoch_control_validation, strict=True)
    ))
    aggregate_target_rate = float(np.sum(chandra_control * target_width))
    aggregate_relative_error = abs(aggregate_target_rate - weighted_native_rate) / max(
        abs(weighted_native_rate), 1.0e-30
    )
    if aggregate_relative_error > 1.0e-10:
        raise RuntimeError(
            "Exposure-weighted center Chandra-RMF control is not count-rate conserving: "
            f"{aggregate_relative_error:.3e}"
        )

    precursor.make_squde_figure(region, squde_data, squde_model, output_dir / "v40_center_squde_data_model_1col.png", 3.35)
    precursor.make_squde_figure(region, squde_data, squde_model, output_dir / "v40_center_squde_data_model_2col.png", 7.00)
    precursor.make_region_resolution_control_figure(
        region,
        squde_data,
        squde_model,
        chandra_control,
        output_dir / "v40_center_squde_mock_chandra_rmf_1col.png",
        3.35,
        control_label="8-epoch Chandra RMF control (SQUDE ARF)",
    )
    precursor.make_region_resolution_control_figure(
        region,
        squde_data,
        squde_model,
        chandra_control,
        output_dir / "v40_center_squde_mock_chandra_rmf_2col.png",
        7.00,
        control_label="8-epoch Chandra RMF control (SQUDE ARF)",
    )

    min50_ylim_top = core.resolution_control_ylim_top(
        squde_data, squde_model, chandra_control
    )
    ungroup(10)
    notice_id(10, 0.3, 4.0)
    group_counts(10, 30)
    grouping_min30 = np.asarray(get_data(10).grouping, dtype=np.int64).copy()
    squde_data_min30 = copy.deepcopy(get_data_plot(10))
    squde_model_min30, squde_components_min30 = fold_components(
        10, average_source, average_components
    )
    squde_closure_min30 = core.validate_component_closure(
        squde_model_min30, squde_components_min30
    )
    explicit_squde_fold_min30 = precursor.validate_explicit_response_fold(
        average_source,
        squde_model_min30,
        precursor.SQUDE_ARF,
        precursor.SQUDE_RMF,
    )
    target_lo_min30 = np.asarray(squde_model_min30.xlo, dtype=float)
    target_hi_min30 = np.asarray(squde_model_min30.xhi, dtype=float)
    explicit_average_min30 = core.rebin_density_conserving_rate(
        squde_kernel.channel_lo,
        squde_kernel.channel_hi,
        core.fold_source_with_kernel(
            average_source, squde_kernel, precursor.SQUDE_ARF
        ),
        target_lo_min30,
        target_hi_min30,
    )
    weighted_epoch_squde_min30 = np.zeros_like(explicit_average_min30)
    for dataset_id, weight in enumerate(weights, 1):
        epoch_density = core.fold_source_with_kernel(
            sources[dataset_id], squde_kernel, precursor.SQUDE_ARF
        )
        weighted_epoch_squde_min30 += weight * core.rebin_density_conserving_rate(
            squde_kernel.channel_lo,
            squde_kernel.channel_hi,
            epoch_density,
            target_lo_min30,
            target_hi_min30,
        )
    linearity_relative_min30 = np.abs(
        weighted_epoch_squde_min30 - explicit_average_min30
    ) / np.maximum(np.abs(explicit_average_min30), 1.0e-30)
    squde_linearity_min30 = {
        "weight_sum": float(np.sum(weights)),
        "median_relative_error": float(np.median(linearity_relative_min30)),
        "p95_relative_error": float(np.percentile(linearity_relative_min30, 95)),
        "max_relative_error": float(np.max(linearity_relative_min30)),
    }
    if squde_linearity_min30["p95_relative_error"] > 1.0e-12:
        raise RuntimeError(
            "Center min30 SQUDE source linearity failed: "
            f"p95={squde_linearity_min30['p95_relative_error']:.3e}"
        )

    epoch_controls_min30: list[np.ndarray] = []
    epoch_control_validation_min30: list[dict[str, Any]] = []
    for dataset_id, epoch in enumerate(epochs, 1):
        epoch_rmf = epoch["pha"].with_suffix(".rmf")
        control, validation = precursor.build_chandra_rmf_control(
            sources[dataset_id],
            epoch_rmf,
            squde_model_min30,
        )
        epoch_controls_min30.append(control)
        epoch_control_validation_min30.append({
            "dataset_id": dataset_id,
            "exposure_weight": float(weights[dataset_id - 1]),
            "rmf": precursor.provenance_path(epoch_rmf),
            "rmf_sha256": sha256(epoch_rmf),
            **validation,
        })
    chandra_control_min30 = np.sum(
        np.asarray(weights)[:, None] * np.asarray(epoch_controls_min30),
        axis=0,
    )
    target_width_min30 = target_hi_min30 - target_lo_min30
    weighted_native_rate_min30 = float(sum(
        weight * validation["native_overlap_rate_counts_s"]
        for weight, validation in zip(
            weights, epoch_control_validation_min30, strict=True
        )
    ))
    aggregate_target_rate_min30 = float(np.sum(
        chandra_control_min30 * target_width_min30
    ))
    aggregate_relative_error_min30 = abs(
        aggregate_target_rate_min30 - weighted_native_rate_min30
    ) / max(abs(weighted_native_rate_min30), 1.0e-30)
    if aggregate_relative_error_min30 > 1.0e-10:
        raise RuntimeError(
            "Center min30 Chandra-RMF control is not count-rate conserving: "
            f"{aggregate_relative_error_min30:.3e}"
        )
    shared_ylim_top = max(
        min50_ylim_top,
        core.resolution_control_ylim_top(
            squde_data_min30, squde_model_min30, chandra_control_min30
        ),
    )
    precursor.make_squde_figure(
        region,
        squde_data_min30,
        squde_model_min30,
        output_dir / "v40_center_squde_data_model_min30_1col.png",
        3.35,
    )
    precursor.make_squde_figure(
        region,
        squde_data_min30,
        squde_model_min30,
        output_dir / "v40_center_squde_data_model_min30_2col.png",
        7.00,
    )
    precursor.make_region_resolution_control_figure(
        region,
        squde_data_min30,
        squde_model_min30,
        chandra_control_min30,
        output_dir / "v40_center_squde_mock_chandra_rmf_min30_1col.png",
        3.35,
        control_label="8-epoch Chandra RMF control (SQUDE ARF)",
        ylim_top=shared_ylim_top,
    )
    slider_geometry_min30 = precursor.make_region_resolution_control_figure(
        region,
        squde_data_min30,
        squde_model_min30,
        chandra_control_min30,
        output_dir / "v40_center_squde_mock_chandra_rmf_min30_2col.png",
        7.00,
        control_label="8-epoch Chandra RMF control (SQUDE ARF)",
        ylim_top=shared_ylim_top,
        slider_canvas=True,
    )
    slider_geometry_min50 = precursor.make_region_resolution_control_figure(
        region,
        squde_data,
        squde_model,
        chandra_control,
        output_dir / "v40_center_squde_mock_chandra_rmf_min50_slider_2col.png",
        7.00,
        control_label="8-epoch Chandra RMF control (SQUDE ARF)",
        ylim_top=shared_ylim_top,
        slider_canvas=True,
    )
    if slider_geometry_min30 != slider_geometry_min50:
        raise RuntimeError("v40_center min30/min50 slider geometry differs")
    slider_geometry = slider_geometry_min30
    slider_pair_id = "v40_center_min30_vs_min50_control_v1"

    thermal_flux = float(calc_energy_flux(0.3, 4.0, model=shared_components["cool_vapec"] + shared_components["hot_vapec"]))
    cx_flux = float(calc_energy_flux(0.3, 4.0, model=shared_components["cx_o"] + shared_components["cx_mg"]))
    xrb_flux = float(calc_energy_flux(0.3, 4.0, model=average_xrb))
    total_flux = thermal_flux + cx_flux + xrb_flux
    count_rates = [epoch["count_rate"] for epoch in epochs]
    parameter_values = fit_parameter_dict(fit_result)

    epoch_records = []
    for epoch in epochs:
        pha = epoch["pha"]
        dataset_id = epoch["dataset_id"]
        arf = pha.with_suffix(".arf")
        rmf = pha.with_suffix(".rmf")
        epoch_records.append({
            "dataset_id": dataset_id,
            "pha": precursor.provenance_path(pha),
            "arf": precursor.provenance_path(arf),
            "rmf": precursor.provenance_path(rmf),
            "pha_sha256": sha256(pha), "arf_sha256": sha256(arf), "rmf_sha256": sha256(rmf),
            "exposure_s": epoch["exposure_s"], "exposure_weight": float(weights[dataset_id - 1]),
            "raw_counts": epoch["raw_counts"], "count_rate": epoch["count_rate"],
            "xrb_nH": float(get_model_component(f"absxrb{dataset_id}").nH.val),
            "xrb_ampl": float(get_model_component(f"xrb{dataset_id}").ampl.val),
            "explicit_rmf_vs_sherpa": epoch_explicit_validation[dataset_id - 1],
        })

    variant_min50 = core.display_grouping_variant_payload(
        squde_data,
        squde_model,
        chandra_control,
        grouping_min50,
        50,
        output_dir,
        {
            "mock_figure_1col": "v40_center_squde_data_model_1col.png",
            "mock_figure_2col": "v40_center_squde_data_model_2col.png",
            "control_figure_1col": "v40_center_squde_mock_chandra_rmf_1col.png",
            "control_figure_2col": "v40_center_squde_mock_chandra_rmf_2col.png",
            "slider_control_2col": "v40_center_squde_mock_chandra_rmf_min50_slider_2col.png",
        },
        slider_geometry,
        slider_pair_id,
    )
    variant_min50["component_closure"] = squde_closure
    variant_min50["explicit_response_vs_sherpa"] = explicit_squde_fold
    variant_min50["source_linearity"] = squde_linearity
    variant_min50["control_rebin_conservation"] = {
        "weighted_native_overlap_rate_counts_s": weighted_native_rate,
        "aggregate_target_rate_counts_s": aggregate_target_rate,
        "aggregate_relative_error": aggregate_relative_error,
        "epochs": epoch_control_validation,
    }
    variant_min30 = core.display_grouping_variant_payload(
        squde_data_min30,
        squde_model_min30,
        chandra_control_min30,
        grouping_min30,
        30,
        output_dir,
        {
            "mock_figure_1col": "v40_center_squde_data_model_min30_1col.png",
            "mock_figure_2col": "v40_center_squde_data_model_min30_2col.png",
            "control_figure_1col": "v40_center_squde_mock_chandra_rmf_min30_1col.png",
            "slider_control_2col": "v40_center_squde_mock_chandra_rmf_min30_2col.png",
        },
        slider_geometry,
        slider_pair_id,
    )
    variant_min30["same_raw_pha"] = squde_pha.name
    variant_min30["component_closure"] = squde_closure_min30
    variant_min30["explicit_response_vs_sherpa"] = explicit_squde_fold_min30
    variant_min30["source_linearity"] = squde_linearity_min30
    variant_min30["control_rebin_conservation"] = {
        "weighted_native_overlap_rate_counts_s": weighted_native_rate_min30,
        "aggregate_target_rate_counts_s": aggregate_target_rate_min30,
        "aggregate_relative_error": aggregate_relative_error_min30,
        "epochs": epoch_control_validation_min30,
    }

    return {
        "key": region.key,
        "title": region.title,
        "geometry": region.geometry,
        "chandra": {
            "fit_mode": "joint_8_epoch_shared_diffuse_epoch_variable_xrb",
            "joint_statval": float(fit_result.statval),
            "joint_dof": int(fit_result.dof),
            "joint_rstat": float(fit_result.rstat),
            "joint_qval": float(fit_result.qval),
            "numpoints": int(fit_result.numpoints),
            "shared_xrb_gamma": float(get_model_component("xrb1").gamma.val),
            "parameter_values": parameter_values,
            "epochs": epoch_records,
            "count_rate_min": min(count_rates), "count_rate_max": max(count_rates),
            "count_rate_ratio": max(count_rates) / min(count_rates),
            "combined_display_pha": precursor.provenance_path(combined_pha),
            "combined_display_pha_sha256": sha256(combined_pha),
            "display_bins": int(len(combined_data.x)),
            "thermal_fraction_0p3_4": thermal_flux / total_flux,
            "xrb_fraction_0p3_4": xrb_flux / total_flux,
            "fraction_estimator": "absorbed source-plane energy-flux fraction from calc_energy_flux over 0.3-4.0 keV; not photon-number flux or detector-count fraction",
            "formal_null_rejected": True,
            "tutorial_only": True,
            "not_final_physical_fit": True,
            "physical_inference_scope": "Improved response-forward tutorial baseline only; not diffuse-only physical inference.",
            "soft_gaussian_semantics": "Phenomenological CX line diagnostics/proxies only; not a complete physical multi-line charge-exchange model.",
            "cx_complete_physical_model": False,
            "component_closure": combined_closure,
            "background_status": "No BACKFILE in legacy PHAs; central background is expected to be weak but remains an explicit limitation.",
            "pileup_status": "Broad-aperture point-source pile-up is not separately modeled; point-source-excised re-extraction remains preferred for diffuse inference.",
            "literature_basis": [
                "https://doi.org/10.1111/j.1365-2966.2008.13128.x",
                "https://doi.org/10.1093/mnrasl/slt145",
                "https://doi.org/10.3847/1538-4357/abc010",
                "https://doi.org/10.3847/0004-637X/829/1/28",
                "https://doi.org/10.1086/511174",
                "https://doi.org/10.1111/j.1365-2966.2011.18466.x",
                "https://doi.org/10.1088/0004-637X/794/1/61",
                "https://doi.org/10.1051/0004-6361/202040209",
            ],
        },
        "squde": {
            "source_semantics": "exposure-weighted average of the eight jointly fitted epoch sources",
            "pha": squde_pha.name,
            "raw_realization": {
                "counts_total": raw_squde_counts,
                "channel_sha256_le_f8": core.array_sha256(raw_channels),
                "counts_sha256_le_f8": core.array_sha256(raw_counts_vector),
                "whole_file_sha256_contract": "not asserted because equivalent FITS-header serialization can differ across output roots",
            },
            "arf": precursor.provenance_path(precursor.SQUDE_ARF), "rmf": precursor.provenance_path(precursor.SQUDE_RMF),
            "arf_sha256": sha256(precursor.SQUDE_ARF), "rmf_sha256": sha256(precursor.SQUDE_RMF),
            "exposure_s": precursor.EXPOSURE_S, "seed": region.seed,
            "display_bins": int(len(squde_data.x)),
            "component_closure": squde_closure,
            "explicit_response_vs_sherpa": explicit_squde_fold,
            "source_linearity": squde_linearity,
            "slider_pair_shared_ylim_top_counts_s_keV": shared_ylim_top,
            "grouping_variants": {
                "min50": variant_min50,
                "min30": variant_min30,
            },
            "resolution_control": {
                "operator": "F_i(E) -> fixed SQUDE ARF -> epoch-specific Chandra RMF_i -> rate-conserving rebin; then exposure-weighted sum over eight epochs",
                "source_semantics": "eight jointly fitted epoch sources, each paired with its own Chandra RMF before exposure weighting",
                "rate_unit": "counts s^-1 keV^-1",
                "display_energy_keV": [0.5, 1.7],
                "weight_basis": "PHA EXPOSURE / sum(PHA EXPOSURE)",
                "weight_sum": float(np.sum(weights)),
                "total_chandra_exposure_s": exposure_sum_s,
                "combined_display_exposure_s": combined_exposure_s,
                "exposure_sum_relative_error": exposure_sum_relative_error,
                "fixed_arf": precursor.provenance_path(precursor.SQUDE_ARF),
                "fixed_arf_sha256": sha256(precursor.SQUDE_ARF),
                "aggregation": "exposure-weighted sum of eight epoch-specific source-times-RMF controls",
                **precursor.resolution_control_payload(squde_data, squde_model, chandra_control),
                "weighted_native_overlap_rate_counts_s": weighted_native_rate,
                "aggregate_target_rate_counts_s": aggregate_target_rate,
                "aggregate_relative_error": aggregate_relative_error,
                "epochs": epoch_control_validation,
                "figure_1col": "v40_center_squde_mock_chandra_rmf_1col.png",
                "figure_2col": "v40_center_squde_mock_chandra_rmf_2col.png",
                "slider_figure_2col": "v40_center_squde_mock_chandra_rmf_min50_slider_2col.png",
                "slider_shared_ylim_top_counts_s_keV": shared_ylim_top,
            },
        },
    }


def main(output_dir: Path | None = None) -> None:
    core.configure_plot_style()
    out = OUT_DIR if output_dir is None else output_dir.expanduser().resolve()
    out.mkdir(parents=True, exist_ok=True)
    fit_rows = precursor.load_fit_rows()
    records = []
    for region in precursor.REGIONS:
        if region.key == "v40_center":
            records.append(process_center(region, out))
        else:
            records.append(precursor.process_region(region, fit_rows, out))
    manifest = {
        "version": "joint-center-v6-complete-grouping-contract",
        "model_semantics": "North/South historical regional baselines; Center eight-epoch shared diffuse plus epoch-variable absorbed XRB/ULX continua, one Fe Kalpha 6.40-keV fluorescence line, and two phenomenological soft CX line proxies.",
        "resolution_control_semantics": "Each regional high-resolution figure is separate from Chandra data and shows SQUDE mock data, the SQUDE-RMF model, and a fixed-SQUDE-ARF Chandra-RMF redistribution control. Center combines eight epoch-specific source-times-RMF controls with exposure weights.",
        "display_grouping_semantics": "min50 and min30 reuse the exact same ungrouped seeded PHA, source models, ARF/RMF and control operators. Only Sherpa display grouping changes from group_counts(50) to group_counts(30); slider assets share axis limits.",
        "soft_gaussian_semantics": "The fixed 0.777/1.234-keV Gaussians are diagnostic phenomenological line proxies, not a complete physical multi-line charge-exchange model.",
        "cx_complete_physical_model": False,
        "formal_null_rejected": True,
        "tutorial_only": True,
        "not_final_physical_fit": True,
        "regions": records,
        "squde_arf_sha256": sha256(precursor.SQUDE_ARF),
        "squde_rmf_sha256": sha256(precursor.SQUDE_RMF),
    }
    path = out / "m82_three_region_joint_center_manifest.json"
    path.write_text(json.dumps(manifest, indent=2, ensure_ascii=False) + "\n", encoding="utf-8")
    print(json.dumps(manifest, indent=2, ensure_ascii=False))
    print(f"Wrote joint-center three-region products to: {out}")


if __name__ == "__main__":
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--output-dir", type=Path)
    main(parser.parse_args().output_dir)