#!/usr/bin/env python
"""Build Chandra and SQUDE tutorial figures for three M82 vertical regions.

Each region produces three scientifically separate figure roles:
1. Real Chandra ACIS data + Chandra-response-folded total and model components.
2. A deterministic 50 ks SQUDE fake_pha + SQUDE-response-folded total.
3. The same SQUDE mock and model plus a fixed-SQUDE-ARF Chandra-RMF
   redistribution control.

North/south use the accepted historical region fits. The pathological historical
center solution is not used; center is refit under an explicit APEC-dominant
tutorial constraint (cool kT <= 2 keV, power-law amplitude <= 1e-3).
"""
from __future__ import annotations

import argparse
import copy
import csv
import hashlib
import json
import os
import sys
from dataclasses import dataclass
from pathlib import Path
from typing import Any

os.environ.setdefault("MPLCONFIGDIR", "/tmp/m82-three-region-mpl")
os.environ.setdefault("XDG_CACHE_HOME", "/tmp/m82-three-region-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,
    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

SPEC_ROOT = ROOT / "data" / "spectra" / "m82_vertical_40x80_80x80_early100ks"
FIT_SUMMARY = ROOT / "data" / "fits" / "m82_vertical_40x80_80x80_early100ks" / "m82_vertical_40x80_80x80_flux_summary.csv"
SQUDE_ARF = ROOT / "SQUDE_response" / "SQUDE_rsp_v1.1.1.arf"
SQUDE_RMF = ROOT / "SQUDE_response" / "SQUDE_rsp_v1.1.2.rmf"
DEFAULT_OUT = ROOT / "data" / "fits" / "m82_three_region_mock_tutorial_precursor_static_center"
EXPOSURE_S = 50000.0
CENTER_PL_MAX = 1.0e-3
CENTER_KT_MAX = 2.0
PARAM_NAMES = [
    "absth_nH", "vcool_kT", "vcool_O", "vcool_Ne", "vcool_Mg",
    "vcool_Si", "vcool_S", "vcool_Fe", "vcool_norm", "vhot_norm",
    "gce_o_norm", "gce_mg_norm", "abspl_nH", "pl_ampl",
]


@dataclass(frozen=True)
class RegionSpec:
    key: str
    title: str
    geometry: str
    seed: int
    fit_mode: str


REGIONS = [
    RegionSpec("v80_north", "North wind", '80"×80"', 2026071901, "historical_region_fit"),
    RegionSpec("v40_center", "Center starburst", '40"×80"', 2026071902, "constrained_apec_dominant_refit"),
    RegionSpec("v80_south", "South wind", '80"×80"', 2026071903, "historical_region_fit"),
]


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 array_sha256(values: np.ndarray) -> str:
    return hashlib.sha256(np.asarray(values, dtype="<f8").tobytes()).hexdigest()


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


def load_fit_rows() -> dict[str, dict[str, str]]:
    with FIT_SUMMARY.open("r", encoding="utf-8", newline="") as handle:
        return {row["bin"]: row for row in csv.DictReader(handle)}


def model_components() -> dict[str, Any]:
    gal = get_model_component("gal")
    absth = get_model_component("absth")
    abspl = get_model_component("abspl")
    return {
        "cool_vapec": gal * absth * get_model_component("vcool"),
        "hot_vapec": gal * absth * get_model_component("vhot"),
        "cx_o": gal * absth * get_model_component("gce_o"),
        "cx_mg": gal * absth * get_model_component("gce_mg"),
        "powerlaw": gal * abspl * get_model_component("pl"),
    }


def apply_frozen_params(params: dict[str, float]) -> core.ModelBundle:
    bundle = core.setup_literature_model(1)
    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")
    abspl = get_model_component("abspl")
    pl = get_model_component("pl")

    assignments = [
        (absth.nH, params["absth_nH"]),
        (vcool.kT, params["vcool_kT"]),
        (vcool.norm, params["vcool_norm"]),
        (vhot.norm, params["vhot_norm"]),
        (gce_o.norm, params["gce_o_norm"]),
        (gce_mg.norm, params["gce_mg_norm"]),
        (abspl.nH, params["abspl_nH"]),
        (pl.ampl, params["pl_ampl"]),
    ]
    for parameter, value in assignments:
        parameter.val = value
        freeze(parameter)
    for element in ["O", "Ne", "Mg", "Si", "S", "Fe"]:
        parameter = getattr(vcool, element)
        parameter.val = params[f"vcool_{element}"]
        freeze(parameter)
    return bundle


def setup_region_model(region: RegionSpec, fit_rows: dict[str, dict[str, str]]) -> tuple[core.ModelBundle, dict[str, Any]]:
    if region.fit_mode == "historical_region_fit":
        row = fit_rows[region.key]
        params = {name: float(row[name]) for name in PARAM_NAMES}
        bundle = apply_frozen_params(params)
        metadata = {
            "fit_mode": region.fit_mode,
            "statval": float(row["statval"]),
            "dof": int(row["dof"]),
            "rstat": float(row["rstat"]),
            "parameter_names": PARAM_NAMES,
            "parameter_values": [params[name] for name in PARAM_NAMES],
            "constraint": None,
        }
        return bundle, metadata

    bundle = core.setup_literature_model(1)
    vcool = get_model_component("vcool")
    pl = get_model_component("pl")
    vcool.kT.max = CENTER_KT_MAX
    pl.ampl = 2.0e-4
    pl.ampl.max = CENTER_PL_MAX
    fit(1)
    result = get_fit_results()
    metadata = {
        "fit_mode": region.fit_mode,
        "statval": float(result.statval),
        "dof": int(result.dof),
        "rstat": float(result.rstat),
        "parameter_names": list(result.parnames),
        "parameter_values": [float(value) for value in result.parvals],
        "constraint": {
            "vcool_kT_max_keV": CENTER_KT_MAX,
            "powerlaw_ampl_max": CENTER_PL_MAX,
            "purpose": "APEC-dominated tutorial baseline; not publication-grade inference",
        },
    }
    return bundle, metadata


def flux_fractions() -> dict[str, float]:
    fluxes = {name: float(calc_energy_flux(0.3, 4.0, model=expr)) for name, expr in model_components().items()}
    total = sum(fluxes.values())
    output = {f"{name}_flux_0p3_4": value for name, value in fluxes.items()}
    output["thermal_fraction"] = (fluxes["cool_vapec"] + fluxes["hot_vapec"]) / total
    output["powerlaw_fraction"] = fluxes["powerlaw"] / total
    return output


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


def make_chandra_figure(
    region: RegionSpec,
    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 = core.COMPONENT_STYLE[name]
        ax.plot(data.x, values, color=color, ls=linestyle, lw=0.8 if width < 4 else 1.05, label=label)
    ax.plot(data.x, model_y, color="#D55E00", lw=1.35, label="Chandra-RMF-folded total model")
    ax.errorbar(data.x, data.y, yerr=data.yerr, fmt=".", ms=2.3, color="0.10", label="Chandra ACIS data")
    ax.set_yscale("log")
    ax.set_xlim(0.3, 4.0)
    ax.set_ylim(ymin, ymax)
    ax.set_xlabel("Energy (keV)")
    ax.set_ylabel(r"Counts s$^{-1}$ keV$^{-1}$")
    ax.set_title(f"M82 {region.title} {region.geometry} Chandra ACIS ({len(data.x)} bins)", fontsize=8.0 if width < 4 else 10.5, pad=58 if width < 4 else 6)
    if width < 4:
        ax.legend(loc="lower center", bbox_to_anchor=(0.5, 1.01), fontsize=4.5, ncol=2, frameon=True, columnspacing=0.8, handlelength=2.4)
    else:
        ax.legend(loc="upper right", fontsize=7.0, ncol=2, frameon=True)
    ax.grid(alpha=0.20)
    fig.savefig(path, dpi=320, bbox_inches="tight")
    plt.close(fig)


def make_squde_figure(region: RegionSpec, data: Any, model: Any, path: Path, width: float) -> None:
    model_y = core.model_on_data_grid(model, data)
    mask = (data.x >= 0.50) & (data.x <= 1.70)
    fig, ax = plt.subplots(figsize=(width, width * 0.62))
    ax.plot(data.x[mask], model_y[mask], color="#D55E00", lw=1.25, label="SQUDE-RMF-folded model")
    ax.errorbar(data.x[mask], data.y[mask], yerr=data.yerr[mask], fmt=".", ms=1.8 if width < 4 else 2.6, color="0.10", label="SQUDE mock data (50 ks)")
    ymax = float(np.max(np.maximum(data.y[mask] + data.yerr[mask], model_y[mask])))
    ax.set_xlim(0.50, 1.70)
    ax.set_ylim(0.0, 1.10 * ymax)
    ax.set_xlabel("Energy (keV)")
    ax.set_ylabel(r"Counts s$^{-1}$ keV$^{-1}$")
    ax.set_title(f"M82 {region.title} {region.geometry}: 50 ks SQUDE mock", fontsize=8.0 if width < 4 else 10.5)
    ax.legend(loc="upper right", fontsize=6.0 if width < 4 else 8.0, frameon=True)
    ax.grid(alpha=0.23)
    fig.savefig(path, dpi=320, bbox_inches="tight")
    plt.close(fig)


def validate_explicit_response_fold(
    source_model: Any,
    native_model: Any,
    chandra_arf: Path,
    chandra_rmf: Path,
) -> dict[str, float]:
    """Require the explicit regional RMF operator to reproduce Sherpa folding."""
    kernel = core.read_rmf_kernel(chandra_rmf)
    explicit_density = core.fold_source_with_kernel(source_model, kernel, chandra_arf)
    explicit = core.rebin_density_conserving_rate(
        kernel.channel_lo,
        kernel.channel_hi,
        explicit_density,
        np.asarray(native_model.xlo, dtype=float),
        np.asarray(native_model.xhi, dtype=float),
    )
    native = np.asarray(native_model.y, dtype=float)
    relative = np.abs(explicit - native) / np.maximum(np.abs(native), 1.0e-30)
    metrics = {
        "median_relative_error": float(np.median(relative)),
        "p95_relative_error": float(np.percentile(relative, 95)),
        "max_relative_error": float(np.max(relative)),
    }
    if metrics["p95_relative_error"] > 1.0e-5:
        raise RuntimeError(
            f"Explicit regional Chandra folding failed for {chandra_rmf}: "
            f"p95={metrics['p95_relative_error']:.3e}"
        )
    return metrics


def build_chandra_rmf_control(
    source_model: Any,
    chandra_rmf: Path,
    target_model: Any,
) -> tuple[np.ndarray, dict[str, float]]:
    """Fold one source through fixed SQUDE ARF + Chandra RMF onto target bins."""
    kernel = core.read_rmf_kernel(chandra_rmf)
    native_density = core.fold_source_with_kernel(source_model, kernel, SQUDE_ARF)
    target_lo = np.asarray(target_model.xlo, dtype=float)
    target_hi = np.asarray(target_model.xhi, dtype=float)
    control = core.rebin_density_conserving_rate(
        kernel.channel_lo,
        kernel.channel_hi,
        native_density,
        target_lo,
        target_hi,
    )
    lower = float(np.min(target_lo))
    upper = float(np.max(target_hi))
    native_overlap = np.clip(
        np.minimum(kernel.channel_hi, upper) - np.maximum(kernel.channel_lo, lower),
        0.0,
        None,
    )
    native_rate = float(np.sum(native_density * native_overlap))
    target_rate = float(np.sum(control * (target_hi - target_lo)))
    relative_error = abs(target_rate - native_rate) / max(abs(native_rate), 1.0e-30)
    if not np.all(np.isfinite(control)) or float(np.min(control)) < -1.0e-12:
        raise RuntimeError(f"Non-finite or negative Chandra-RMF control for {chandra_rmf}")
    if relative_error > 1.0e-10:
        raise RuntimeError(
            f"Chandra-RMF control rebin is not count-rate conserving: {relative_error:.3e}"
        )
    return control, {
        "native_overlap_rate_counts_s": native_rate,
        "target_rate_counts_s": target_rate,
        "relative_error": relative_error,
    }


def resolution_control_payload(
    data: Any,
    squde_model: Any,
    chandra_control: np.ndarray,
) -> dict[str, Any]:
    """Serialize the exact model curves plotted on the SQUDE data grid."""
    squde_display = core.model_on_data_grid(squde_model, data)
    chandra_display = core.rebin_density_conserving_rate(
        np.asarray(squde_model.xlo, dtype=float),
        np.asarray(squde_model.xhi, dtype=float),
        np.asarray(chandra_control, dtype=float),
        np.asarray(data.xlo, dtype=float),
        np.asarray(data.xhi, dtype=float),
    )
    return {
        "model_grid_chandra_curve_sha256_float64_le": array_sha256(chandra_control),
        "display_grid": {
            "energy_lo_keV": [float(value) for value in np.asarray(data.xlo)],
            "energy_hi_keV": [float(value) for value in np.asarray(data.xhi)],
            "squde_rmf_counts_s_keV": [float(value) for value in squde_display],
            "chandra_rmf_squde_arf_counts_s_keV": [float(value) for value in chandra_display],
            "squde_curve_sha256_float64_le": array_sha256(squde_display),
            "chandra_curve_sha256_float64_le": array_sha256(chandra_display),
        },
    }


def make_region_resolution_control_figure(
    region: RegionSpec,
    data: Any,
    squde_model: Any,
    chandra_control: np.ndarray,
    path: Path,
    width: float,
    control_label: str = "Chandra RMF model (SQUDE ARF)",
    ylim_top: float | None = None,
    slider_canvas: bool = False,
) -> dict[str, Any]:
    return core.make_resolution_control_figure(
        data,
        squde_model,
        chandra_control,
        path,
        width,
        title=f"M82 {region.title} {region.geometry}: SQUDE versus Chandra RMF",
        control_label=control_label,
        ylim_top=ylim_top,
        slider_canvas=slider_canvas,
    )


def process_region(region: RegionSpec, fit_rows: dict[str, dict[str, str]], out_dir: Path) -> dict[str, Any]:
    clean()
    pha = SPEC_ROOT / region.key / f"{region.key}_early100ks_combined_src.pi"
    arf = pha.with_suffix(".arf")
    rmf = pha.with_suffix(".rmf")
    load_pha(1, str(pha))
    set_analysis(1, "energy")
    notice_id(1, 0.3, 7.0)
    group_counts(1, 20)
    set_stat("chi2gehrels")
    set_method("levmar")
    bundle, fit_metadata = setup_region_model(region, fit_rows)
    fractions = flux_fractions()

    chandra_data = copy.deepcopy(get_data_plot(1))
    chandra_model = copy.deepcopy(get_model_plot(1))
    chandra_components = core.response_folded_components(bundle, 1)
    chandra_model = copy.deepcopy(get_model_plot(1))
    chandra_closure = core.validate_component_closure(chandra_model, chandra_components)
    explicit_chandra_fold = validate_explicit_response_fold(
        bundle.total,
        chandra_model,
        arf,
        rmf,
    )

    make_chandra_figure(region, chandra_data, chandra_model, chandra_components, out_dir / f"{region.key}_chandra_data_model_components_1col.png", 3.35)
    make_chandra_figure(region, chandra_data, chandra_model, chandra_components, out_dir / f"{region.key}_chandra_data_model_components_2col.png", 7.00)

    set_rng(np.random.default_rng(region.seed))
    set_source(2, bundle.total)
    squde_pha = out_dir / f"{region.key}_squde_50ks.pha"
    fake_pha(2, arf=str(SQUDE_ARF), rmf=str(SQUDE_RMF), exposure=EXPOSURE_S, grouped=False, method=poisson_noise)
    save_pha(2, str(squde_pha), clobber=True)
    raw_channels = np.asarray(get_data(2).channel, dtype=float).copy()
    raw_counts_vector = np.asarray(get_data(2).counts, dtype=float).copy()
    raw_squde_counts = int(np.sum(raw_counts_vector))
    set_analysis(2, "energy")
    notice_id(2, 0.3, 4.0)
    group_counts(2, 50)
    grouping_min50 = np.asarray(get_data(2).grouping, dtype=np.int64).copy()
    squde_data = copy.deepcopy(get_data_plot(2))
    squde_components = core.response_folded_components(bundle, 2)
    squde_model = copy.deepcopy(get_model_plot(2))
    squde_closure = core.validate_component_closure(squde_model, squde_components)
    explicit_squde_fold = validate_explicit_response_fold(
        bundle.total,
        squde_model,
        SQUDE_ARF,
        SQUDE_RMF,
    )
    chandra_control, control_validation = build_chandra_rmf_control(
        bundle.total,
        rmf,
        squde_model,
    )

    make_squde_figure(region, squde_data, squde_model, out_dir / f"{region.key}_squde_data_model_1col.png", 3.35)
    make_squde_figure(region, squde_data, squde_model, out_dir / f"{region.key}_squde_data_model_2col.png", 7.00)
    make_region_resolution_control_figure(
        region,
        squde_data,
        squde_model,
        chandra_control,
        out_dir / f"{region.key}_squde_mock_chandra_rmf_1col.png",
        3.35,
    )
    make_region_resolution_control_figure(
        region,
        squde_data,
        squde_model,
        chandra_control,
        out_dir / f"{region.key}_squde_mock_chandra_rmf_2col.png",
        7.00,
    )

    min50_ylim_top = core.resolution_control_ylim_top(
        squde_data, squde_model, chandra_control
    )
    ungroup(2)
    notice_id(2, 0.3, 4.0)
    group_counts(2, 30)
    grouping_min30 = np.asarray(get_data(2).grouping, dtype=np.int64).copy()
    squde_data_min30 = copy.deepcopy(get_data_plot(2))
    squde_components_min30 = core.response_folded_components(bundle, 2)
    squde_model_min30 = copy.deepcopy(get_model_plot(2))
    squde_closure_min30 = core.validate_component_closure(
        squde_model_min30, squde_components_min30
    )
    explicit_squde_fold_min30 = validate_explicit_response_fold(
        bundle.total,
        squde_model_min30,
        SQUDE_ARF,
        SQUDE_RMF,
    )
    chandra_control_min30, control_validation_min30 = build_chandra_rmf_control(
        bundle.total,
        rmf,
        squde_model_min30,
    )
    shared_ylim_top = max(
        min50_ylim_top,
        core.resolution_control_ylim_top(
            squde_data_min30, squde_model_min30, chandra_control_min30
        ),
    )
    make_squde_figure(
        region,
        squde_data_min30,
        squde_model_min30,
        out_dir / f"{region.key}_squde_data_model_min30_1col.png",
        3.35,
    )
    make_squde_figure(
        region,
        squde_data_min30,
        squde_model_min30,
        out_dir / f"{region.key}_squde_data_model_min30_2col.png",
        7.00,
    )
    make_region_resolution_control_figure(
        region,
        squde_data_min30,
        squde_model_min30,
        chandra_control_min30,
        out_dir / f"{region.key}_squde_mock_chandra_rmf_min30_1col.png",
        3.35,
        ylim_top=shared_ylim_top,
    )
    slider_geometry_min30 = make_region_resolution_control_figure(
        region,
        squde_data_min30,
        squde_model_min30,
        chandra_control_min30,
        out_dir / f"{region.key}_squde_mock_chandra_rmf_min30_2col.png",
        7.00,
        ylim_top=shared_ylim_top,
        slider_canvas=True,
    )
    slider_geometry_min50 = make_region_resolution_control_figure(
        region,
        squde_data,
        squde_model,
        chandra_control,
        out_dir / f"{region.key}_squde_mock_chandra_rmf_min50_slider_2col.png",
        7.00,
        ylim_top=shared_ylim_top,
        slider_canvas=True,
    )
    if slider_geometry_min30 != slider_geometry_min50:
        raise RuntimeError(f"{region.key} min30/min50 slider geometry differs")
    slider_geometry = slider_geometry_min30
    slider_pair_id = f"{region.key}_min30_vs_min50_control_v1"
    variant_min50 = core.display_grouping_variant_payload(
        squde_data,
        squde_model,
        chandra_control,
        grouping_min50,
        50,
        out_dir,
        {
            "mock_figure_1col": f"{region.key}_squde_data_model_1col.png",
            "mock_figure_2col": f"{region.key}_squde_data_model_2col.png",
            "control_figure_1col": f"{region.key}_squde_mock_chandra_rmf_1col.png",
            "control_figure_2col": f"{region.key}_squde_mock_chandra_rmf_2col.png",
            "slider_control_2col": f"{region.key}_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["control_rebin_conservation"] = control_validation
    variant_min30 = core.display_grouping_variant_payload(
        squde_data_min30,
        squde_model_min30,
        chandra_control_min30,
        grouping_min30,
        30,
        out_dir,
        {
            "mock_figure_1col": f"{region.key}_squde_data_model_min30_1col.png",
            "mock_figure_2col": f"{region.key}_squde_data_model_min30_2col.png",
            "control_figure_1col": f"{region.key}_squde_mock_chandra_rmf_min30_1col.png",
            "slider_control_2col": f"{region.key}_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["control_rebin_conservation"] = control_validation_min30

    return {
        "key": region.key,
        "title": region.title,
        "geometry": region.geometry,
        "chandra": {
            "pha": provenance_path(pha), "arf": provenance_path(arf), "rmf": provenance_path(rmf),
            "pha_sha256": sha256(pha), "arf_sha256": sha256(arf), "rmf_sha256": sha256(rmf),
            "exposure_s": float(get_data(1).exposure), "display_bins": int(len(chandra_data.x)),
            "fit": fit_metadata, "flux_fractions": fractions,
            "component_closure": chandra_closure,
            "explicit_rmf_vs_sherpa": explicit_chandra_fold,
        },
        "squde": {
            "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": provenance_path(SQUDE_ARF), "rmf": provenance_path(SQUDE_RMF),
            "arf_sha256": sha256(SQUDE_ARF), "rmf_sha256": sha256(SQUDE_RMF),
            "exposure_s": EXPOSURE_S, "seed": region.seed,
            "display_bins": int(len(squde_data.x)), "component_closure": squde_closure,
            "explicit_response_vs_sherpa": explicit_squde_fold,
            "slider_pair_shared_ylim_top_counts_s_keV": shared_ylim_top,
            "grouping_variants": {
                "min50": variant_min50,
                "min30": variant_min30,
            },
            "resolution_control": {
                "operator": "F(E) -> fixed SQUDE ARF -> region combined Chandra RMF -> rate-conserving rebin to SQUDE model bins",
                "source_semantics": "the exact fitted regional source used to generate the SQUDE mock",
                "rate_unit": "counts s^-1 keV^-1",
                "display_energy_keV": [0.5, 1.7],
                "chandra_rmf": provenance_path(rmf),
                "chandra_rmf_sha256": sha256(rmf),
                "fixed_arf": provenance_path(SQUDE_ARF),
                "fixed_arf_sha256": sha256(SQUDE_ARF),
                "aggregation": "single region-specific combined Chandra RMF",
                **resolution_control_payload(squde_data, squde_model, chandra_control),
                "rebin_conservation": control_validation,
                "figure_1col": f"{region.key}_squde_mock_chandra_rmf_1col.png",
                "figure_2col": f"{region.key}_squde_mock_chandra_rmf_2col.png",
                "slider_figure_2col": f"{region.key}_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_dir = DEFAULT_OUT if output_dir is None else output_dir.expanduser().resolve()
    out_dir.mkdir(parents=True, exist_ok=True)
    fit_rows = load_fit_rows()
    manifest = {
        "model": "tbabs_gal * (tbabs_wind * (vapec_cool + vapec_hot + gaussian_0.777 + gaussian_1.234) + tbabs_powerlaw * powerlaw)",
        "regions": [process_region(region, fit_rows, out_dir) for region in REGIONS],
        "squde_arf_sha256": sha256(SQUDE_ARF),
        "squde_rmf_sha256": sha256(SQUDE_RMF),
    }
    (out_dir / "m82_three_region_mock_tutorial_manifest.json").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 three-region tutorial products to: {out_dir}")


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