simulation.py

Go to the documentation of this file.

"""
simulation.py — RadCluster_2_1 simulation orchestrator.
 
Responsibilities
----------------
1. Load inputs and initialise physics modules.
2. Dispatch to the appropriate solver mode and physics option.
3. Write timestamped output directory with provenance, CSV, and plots.
 
Solver modes
------------
full_system   → C++ SUNDIALS CVODE BDF, full system, via cpp_bridge.run_cpp_solver
active_window → C++ SUNDIALS CVODE BDF with two independent sliding windows
                (SIA + VAC) and OpenMP-parallel RHS. Thread count is auto-set
                from N_eq; falls back to single-threaded sliding when OpenMP
                is unavailable or only one thread is selected.
 
Legacy aliases (still accepted): cpp_full → full_system,
                                 sliding_OpenMP → active_window.
 
Physics options (two orthogonal axes)
-------------------------------------
equations ∈ {'discrete', 'bin_moment'} : per-size ODEs vs. bin-moment grouping
cascade   ∈ {'fission',  'fusion'}     : He coupling case + cascade spectrum
 
Combinations map to the canonical physics_option string:
  discrete   + fission → 'full_CD_fission'        (he_mode='case2', Eq. 175)
  discrete   + fusion  → 'full_CD_fusion'         (he_mode='case1', Eq. 174)
  bin_moment + fission → 'bin_moment_CD_fission'  (he_mode='case2')
  bin_moment + fusion  → 'bin_moment_CD_fusion'   (he_mode='case1')
 
Either form is accepted by RadClusterSimulation.  Legacy equations='full_CD'
is silently aliased to 'discrete'.
"""
 
import time as _time
import os
import platform
import re
import signal as _signal
import sys
from pathlib import Path
from datetime import datetime
 
import numpy as np
 
from .input_data      import InputData
from .reaction_rates  import ReactionRates
from .rate_equations  import RateEquations
from .bin_moment_rates import BinMomentRateEquations
from .                import post_process
 
 
BASE_DIR = Path(__file__).parent.parent
 
 
def _remap_bin_moments(mu0_old, mu1_old, bins_old, bins_new):
    """Move (μ₀, μ₁) from one bin partition to another, conserving both sums.
 
    Used when a domain doubling re-derives the geometric bin edges and the
    moments have to be expressed on the new partition.
 
    WHY NOT reconstruct-then-reproject.  The obvious route — expand the moments
    to a per-size c_n with the active shape function, then re-take moments on
    the new bins — is NOT conservative.  ``distribution_from_moments_hat``
    clamps its negative lobe with ``np.maximum(vals, 0.0)``, and that clip adds
    mass: for a bin whose mean sits away from the centre the reconstructed sum
    exceeds the bin's own μ₀, without bound.  Measured on the ⟨100⟩ block of the
    I=40000 run (2026-08-15): Σ(reconstructed) / Σμ₀ = 176.9, which is precisely
    the ×177 jump in N_100 across that run's doubling, and it took δ_FP from
    0.20 to 0.96.  The ½⟨111⟩ and vacancy blocks survive the same round trip
    only because they are broad and well-centred; the ⟨100⟩ population is
    sessile, one-way and sharply peaked, so it sits in the worst case
    permanently.
 
    THIS route never builds a per-size array.  Each old bin is treated as
    uniformly filled and split across the new bins it overlaps, by overlap
    fraction w (Σ_K w = 1 because the new bins tile the old bin's span).  Both
    moments ride the same weights, so Σμ₀ and Σμ₁ are conserved to round-off
    by construction — the clip cannot enter because no reconstruction happens.
 
    The one approximation is the within-bin placement of content: splitting μ₁
    by width smears a peak across the pieces of a straddled bin.  A final clamp
    puts each new bin's mean back inside its own range, mirroring the
    consistency rule in ``post_process._floor_bin_moments`` (an empty bin
    carries no content; μ₁/μ₀ lies in the bin).  Both are second-order next to
    a factor of 177.
 
    Returns (mu0_new, mu1_new), each length len(bins_new).
    """
    mu0_old = np.asarray(mu0_old, dtype=float)
    K_new = len(bins_new)
    mu0_new = np.zeros(K_new)
    mu1_new = np.zeros(K_new)
    has_m1 = mu1_old is not None
    mu1_old = (np.asarray(mu1_old, dtype=float) if has_m1
               else mu0_old * np.array([0.5 * (lo + hi - 1)
                                        for lo, hi in bins_old]))
 
    for k, (lo, hi) in enumerate(bins_old):
        width = float(hi - lo)
        if width <= 0.0 or mu0_old[k] <= 0.0:
            continue
        for K, (Lo, Hi) in enumerate(bins_new):
            ov = min(hi, Hi) - max(lo, Lo)
            if ov <= 0:
                continue
            w = ov / width
            mu0_new[K] += w * mu0_old[k]
            mu1_new[K] += w * mu1_old[k]
 
    # Consistency: empty bin carries no content; mean size lies in the bin.
    for K, (Lo, Hi) in enumerate(bins_new):
        if mu0_new[K] <= 0.0:
            mu0_new[K] = 0.0
            mu1_new[K] = 0.0
            continue
        nb = mu1_new[K] / mu0_new[K]
        if nb < Lo:
            mu1_new[K] = mu0_new[K] * float(Lo)
        elif nb > Hi - 1:
            mu1_new[K] = mu0_new[K] * float(Hi - 1)
    return mu0_new, mu1_new
 
 
class _InterruptDeferred:
    """Defer SIGINT through a critical save block.
 
    First Ctrl+C inside the block prints a notice and continues — the save
    finishes flushing partial output, then KeyboardInterrupt is re-raised
    on exit.  A second Ctrl+C inside the same block restores the default
    handler so the user can force-abort if a save step is hung.
 
    Only used in main-thread Python contexts; signal.signal raises
    ValueError off the main thread, so we silently fall back to a no-op.
    """
    __slots__ = ('_msg', '_caught', '_prev', '_active')
 
    def __init__(self, msg=('Ctrl+C — finishing pending writes before exit. '
                            'Press Ctrl+C again to force-abort.')):
        self._msg     = msg
        self._caught  = False
        self._prev    = None
        self._active  = False
 
    def __enter__(self):
        try:
            self._prev = _signal.getsignal(_signal.SIGINT)
            _signal.signal(_signal.SIGINT, self._handler_handler)
            self._active = True
        except (ValueError, OSError):
            # Off main thread, or signal not available — no-op.
            self._active = False
        return self
 
    def _handler(self, signum, frame):
        if not self._caught:
            self._caught = True
            try:
                print(f"\n*** {self._msg}", flush=True)
            except Exception:
                pass
        else:
            # User insisted — restore previous handler and re-raise now.
            try:
                _signal.signal(_signal.SIGINT, self._prev)
            except Exception:
                pass
            self._active = False
            raise KeyboardInterrupt
 
    def __exit__(self, exc_type, exc, tb):
        if self._active:
            try:
                _signal.signal(_signal.SIGINT, self._prev)
            except Exception:
                pass
            self._active = False
        # If a Ctrl+C arrived during the block and no other exception is
        # already propagating, raise it now so the caller still sees the
        # interrupt — just after the partial output has been written.
        if self._caught and exc_type is None:
            raise KeyboardInterrupt
        return False
 
 
class RadClusterSimulation:
    """
    Orchestrates the RadCluster_2_1 simulation workflow.
 
    Parameters
    ----------
    I              : int, optional  — override max SIA cluster size
    V              : int, optional  — override max vacancy cluster size
    solver_mode    : str, optional  — 'full_system' | 'active_window'
    physics_option : str, optional  — 'full_CD_fission' | 'full_CD_fusion' |
                                      'bin_moment_CD_fission' | 'bin_moment_CD_fusion'.
                                      Mutually exclusive with the new (equations,
                                      cascade) pair below.
    equations      : str, optional  — 'discrete' | 'bin_moment'.  Use together with
                                      `cascade` instead of `physics_option`.
                                      Legacy 'full_CD' is silently aliased to
                                      'discrete'.
    cascade        : str, optional  — 'fission' | 'fusion'.  Use together with
                                      `equations`.
    excel_file     : path-like, optional
    C_floor        : float, optional — concentration floor (default 1e-15).
                                       Any state variable below this is clamped
                                       before evaluating rate terms.
    he_kinetics    : str, optional   — 'dynamic' (default) integrates free He as a
                                       full ODE (Eq. 157); 'quasi_steady_state'
                                       eliminates it via dc_h/dt = 0 (valid because
                                       E_m_h = 0.06 eV gives rapid equilibration).
    i_mobile       : int, optional   — max mobile SIA cluster size (default 1).
                                       1 = mono-defect form; set > 1 to enable 1D
                                       glide and SIA–SIA coalescence.
    v_mobile       : int, optional   — max mobile vacancy cluster size (default 1).
                                       1 = mono-defect form; set > 1 to enable
                                       vacancy–vacancy coalescence.
 
    Notes
    -----
    Two-layer RAG description.  Alongside the legacy RateEquations /
    BinMomentRateEquations objects that drive the C++ solver, every
    RadClusterSimulation also carries the abstract-core + EUROFER-97
    instantiation as ``self.material_rag`` (the reaction admissibility
    graph G_Eur, paper Section 4) and ``self.state_layout`` (its
    contiguous ODE-index map).  This makes the graph-based formulation
    the live structural description of every run.  The RAG build is
    advisory: it is wrapped in try/except so a RAG-build issue can never
    break a run, and the C++ drive path is byte-identical with or
    without it.  ``material_rag`` / ``state_layout`` are ``None`` if the
    build failed.
    """
 
    def __init__(self, I=None, V=None, solver_mode=None,
                 physics_option=None, equations=None, cascade=None,
                 excel_file=None,
                 C_floor=None, he_kinetics=None,
                 i_mobile=None, v_mobile=None,
                 **legacy_kw):
        # Backward compatibility: accept old N/M and n_max_i/m_max_v kwargs
        if I is None and 'N' in legacy_kw:
            I = legacy_kw.pop('N')
        if V is None and 'M' in legacy_kw:
            V = legacy_kw.pop('M')
        if i_mobile is None and 'n_max_i' in legacy_kw:
            i_mobile = legacy_kw.pop('n_max_i')
        if v_mobile is None and 'm_max_v' in legacy_kw:
            v_mobile = legacy_kw.pop('m_max_v')
 
        # Anything still in legacy_kw is a name this constructor does not
        # know.  Silently dropping it is dangerous: on 2026-08-02 a nine-case
        # bin-moment refinement study passed i_discrete=/I_bin=/shape_function=
        # here, had all nine silently ignored, produced nine bit-identical
        # runs, and was very nearly reported as a structural failure of the
        # bin-moment closure.  Fail loudly instead.
        if legacy_kw:
            _BIN_CFG = {'i_discrete', 'v_discrete', 'I_bin', 'V_bin',
                        'shape_function', 'bin_ratio'}
            bad = sorted(legacy_kw)
            msg = (f"RadClusterSimulation() got unexpected keyword argument(s): "
                   f"{', '.join(bad)}")
            if _BIN_CFG.intersection(bad):
                msg += (
                    "\n  Bin-moment configuration is NOT a constructor argument. "
                    "Set it on the workbook dict and rebuild:\n"
                    "      sim.input_data.reactions['i_discrete'] = ...\n"
                    "      sim.input_data._calculate_derived(); sim.rebuild_rates()\n"
                    "  Read the values back from sim.rate_equations (NOT "
                    "sim.reaction_rates) to confirm they took effect.")
            raise TypeError(msg)
 
        # Resolve (equations, cascade) → physics_option.  Either form is
        # accepted, but not both at once.
        if equations is not None or cascade is not None:
            if physics_option is not None:
                raise ValueError(
                    "Pass either physics_option= OR (equations=, cascade=), "
                    "not both.")
            if equations is None or cascade is None:
                raise ValueError(
                    "Both equations= and cascade= must be provided together "
                    "(or use physics_option= instead).")
            from .input_data import make_physics_option
            physics_option = make_physics_option(equations, cascade)
 
        print("Initializing RadCluster_2_1 simulation…")
        kwargs = {}
        if excel_file is not None:
            kwargs['excel_file'] = excel_file
 
        self.input_data = InputData(
            I=I, V=V,
            solver_mode=solver_mode,
            physics_option=physics_option,
            **kwargs
        )
 
        # Inject optional overrides before building rate equations
        if i_mobile is not None:
            self.input_data.diffusion['i_mobile'] = int(i_mobile)
            self.input_data.derived['i_mobile']   = int(i_mobile)
        if v_mobile is not None:
            self.input_data.diffusion['v_mobile'] = int(v_mobile)
            self.input_data.derived['v_mobile']   = int(v_mobile)
        if C_floor is not None:
            self.input_data.reactions['C_floor'] = float(C_floor)
        if he_kinetics is not None:
            validated = str(he_kinetics).lower()
            if validated not in ('dynamic', 'quasi_steady_state'):
                import warnings
                warnings.warn(f"Unknown he_kinetics='{he_kinetics}'. Using 'dynamic'.")
                validated = 'dynamic'
            self.input_data.reactions['he_kinetics'] = validated
 
        self.reaction_rates = ReactionRates(self.input_data)
 
        po = self.input_data.physics_option
        if 'bin_moment' in po:
            self.rate_equations = BinMomentRateEquations(
                self.input_data, self.reaction_rates
            )
        else:
            self.rate_equations = RateEquations(
                self.input_data, self.reaction_rates
            )
 
        # ── Two-layer RAG description (abstract core + EUROFER-97) ───────────
        # ADDITIVE: build the EUROFER-97 reaction admissibility graph and its
        # StateLayout alongside the legacy RateEquations object above.  This
        # does NOT touch the C++ drive path — it only attaches the graph-based
        # structural description of the run to the simulation object.  Wrapped
        # in try/except so a RAG-build issue can never break a run.
        self.material_rag = None
        self.state_layout = None
        try:
            from .materials.eurofer97 import build_eurofer_rag
            self.material_rag, self.state_layout = build_eurofer_rag(
                self.input_data, self.reaction_rates,
                equations=self.input_data.equations,
                cascade=self.input_data.cascade,
            )
            print(f"Material RAG built: {self.material_rag.summary()}")
        except Exception as exc:  # noqa: BLE001 — advisory build, never fatal
            import warnings
            warnings.warn(
                f"EUROFER-97 RAG build failed ({exc!r}); continuing without "
                "the two-layer structural description. The run is unaffected.")
 
        # Always defined so notebook fallbacks like
        # `sim._accumulated_results` don't AttributeError when run_adaptive
        # is interrupted before it can initialize this attribute itself.
        self._accumulated_results = None
 
        print(f"Simulation initialized: solver_mode='{self.input_data.solver_mode}'"
              f"  equations='{self.input_data.equations}'"
              f"  cascade='{self.input_data.cascade}'"
              f"  (physics_option='{po}')")
 
    def rebuild_rates(self):
        """Rebuild reaction rates and rate equations from current input_data.
 
        Call this after modifying input_data parameter dicts and calling
        input_data._calculate_derived() to ensure the solver uses the
        updated rate constants.
        """
        self.reaction_rates = ReactionRates(self.input_data)
        po = self.input_data.physics_option
        if 'bin_moment' in po:
            self.rate_equations = BinMomentRateEquations(
                self.input_data, self.reaction_rates
            )
        else:
            self.rate_equations = RateEquations(
                self.input_data, self.reaction_rates
            )
 
    def _reconstruct_sia_per_size(self, sia_block):
        """Reconstruct the per-size ½⟨111⟩ distribution c_n (n = 1..I) from a
        SIA state block.
 
        In discrete mode the block is already per-size and is returned (clipped)
        as-is.  In bin-moment mode the block is ``i_discrete`` discrete sizes
        followed by ``I_bin`` bins of ``P`` moments each; the binned part is
        expanded through the active intra-bin shape function, exactly as
        cpp_bridge reconstructs the appended ⟨100⟩ block.
        """
        re_obj = self.rate_equations
        I = self.input_data.I
        col = np.maximum(np.asarray(sia_block, dtype=float), 0.0)
        if 'bin_moment' not in self.input_data.physics_option:
            c = np.zeros(I)
            c[:min(I, col.size)] = col[:I]
            return c
        from .bin_moment_rates import reconstruct_distribution
        i_d  = int(getattr(re_obj, 'i_discrete', 0))
        P    = int(getattr(re_obj, 'n_mom', 2))
        bins = getattr(re_obj, 'bins', [])
        I_bin = len(bins)
        c = np.zeros(I)
        c[:i_d] = col[:i_d]
        if I_bin > 0:
            mom = col[i_d:i_d + P * I_bin]
            mu0 = mom[0::P][:I_bin]
            mu1 = mom[1::P][:I_bin] if P >= 2 else None
            mu2 = mom[2::P][:I_bin] if P >= 3 else None
            rec = reconstruct_distribution(
                getattr(re_obj, 'shape_function', 'linear'),
                mu0, mu1, mu2, bins, I)
            c[i_d:] = rec[i_d:]
        return c
 
    def _loop_network_update(self, results, dt):
        """Operator-split advance of the dynamic network density ρ_net between
        integration segments (loop_network_loss.tex §Implementation).
 
        ρ_net is held piecewise-constant *within* a segment (so the Jacobian /
        Woodbury structure is untouched) and advanced *between* segments by an
        explicit Euler step of Eq. (rho_net_balance):
 
            ρ_net ← ρ_net + dt · [ (π/Ω) Σ d_n Λ_n c_n  −  K_rec ρ_net^{3/2} ]
 
        The segment-end SIA/vacancy monomers are frozen into ``ci1_seg`` /
        ``cv1_seg`` (they set the network climb velocity v_net for the next
        segment), and ``rebuild_rates`` refreshes ReactionRates so the P4 sink
        and Λ_n^net pick up the new ρ_net.  No-op unless the channel is on.
 
        The per-size loop distributions that drive the growth term come from the
        post-processed segment and so work in BOTH discrete and bin-moment modes:
          * ½⟨111⟩ — the SIA block of ``results['y']`` (per-size in discrete
            mode; reconstructed from the SIA bin moments via the active shape
            function in bin-moment mode);
          * ⟨100⟩  — ``results['y_sia100']``, which cpp_bridge has already
            reconstructed to per-size [I, n_t] in both modes.  (Note ``results['y']``
            no longer carries the appended ⟨100⟩ block — cpp_bridge splits it
            off — so it must be read from this dedicated key, not the y tail.)
        """
        rr = self.reaction_rates
        if not getattr(rr, 'loop_network_loss', False) or dt <= 0.0:
            return
        re_obj = self.rate_equations
        I = self.input_data.I
        y_end = np.asarray(results['y'][:, -1], dtype=float)
        i_vac = int(getattr(re_obj, 'i_VAC', I))
        ci1 = float(y_end[0]) if y_end.size > 0 else 0.0
        cv1 = float(y_end[i_vac]) if y_end.size > i_vac else 0.0
 
        # Per-size ½⟨111⟩ loop distribution (the SIA block precedes the VAC block).
        c_111 = self._reconstruct_sia_per_size(y_end[:i_vac])
 
        # Per-size ⟨100⟩ loop distribution — already reconstructed by cpp_bridge
        # (per-size [I, n_t] in both discrete and bin-moment modes).
        y100 = results.get('y_sia100')
        c_100 = None
        if y100 is not None and np.asarray(y100).size:
            c_100 = np.asarray(y100, dtype=float)[:, -1]
 
        drho = rr.loop_network_drho_dt(c_111, c_100)
        rho_floor = float(self.input_data.reactions.get('rho_d', 1.0e14))
        # Ceiling guards against ρ_net runaway when K_rec is small/uncalibrated:
        # the gain term grows ∝ ρ_net (Λ ∝ ρ_net), so without the mandatory
        # recovery sink ρ_net would diverge and drive the sinks numerically
        # stiff.  Default ceiling = a physically large network density.
        from .reaction_rates import _num
        rho_max = _num(self.input_data.reactions.get('loop_net_rho_max', 1.0e16),
                       1.0e16)
        rho_new = min(max(float(rr.rho_net) + dt * drho, rho_floor), rho_max)
 
        re = self.input_data.reactions
        re['rho_net'] = rho_new
        re['ci1_seg'] = ci1
        re['cv1_seg'] = cv1
        self.rebuild_rates()
 
    def _resize_domain(self, I_new, V_new):
        """Resize the cluster domain and rebuild all rate objects."""
        self.input_data.reactions['I'] = int(I_new)
        self.input_data.reactions['V'] = int(V_new)
        # InputData caches I, V as properties from reactions dict
        self.rebuild_rates()
 
    def _boundary_fraction_at(self, results, idx):
        """
        Compute tail fractions at a specific time index.
 
        Returns (frac_I, frac_V) — fraction of the total SIA/vacancy content
        in the last 10% of sizes.  Values > ~0.01 indicate boundary pileup.
        """
        y = results.get('y')
        if y is None:
            return 0.0, 0.0
 
        re = self.rate_equations
        is_bin = hasattr(re, 'bins')
        I, V = self.input_data.I, self.input_data.V
        ns = np.arange(1, I + 1, dtype=float)
        ms = np.arange(1, V + 1, dtype=float)
 
        yj = np.maximum(y[:, idx], 0.0)
 
        if is_bin:
            from .bin_moment_rates import reconstruct_distribution
            I_bin = getattr(re, 'I_bin', getattr(re, 'K', 0))
            V_bin = getattr(re, 'V_bin', getattr(re, 'K_v', 0))
            P     = getattr(re, 'n_mom', 2)
            sf    = getattr(re, 'shape_function', 'linear')
            i_d   = getattr(re, 'i_discrete', 0)
            v_d   = getattr(re, 'v_discrete', 0)
            iv    = getattr(re, 'i_VAC', i_d + P * I_bin)
 
            # SIA: discrete + binned reconstruction
            if I_bin > 0:
                mom = yj[i_d:i_d + P * I_bin]
                mu0 = mom[0::P][:I_bin]
                mu1 = mom[1::P][:I_bin] if P >= 2 else None
                mu2 = mom[2::P][:I_bin] if P >= 3 else None
                c_n = reconstruct_distribution(sf, mu0, mu1, mu2, re.bins, I)
                c_n[:i_d] = yj[:i_d]
            else:
                c_n = np.zeros(I)
                c_n[:i_d] = yj[:i_d]
 
            # VAC: discrete + binned reconstruction
            if V_bin > 0:
                vac_start = iv + v_d
                vmom = yj[vac_start:vac_start + P * V_bin]
                vmu0 = vmom[0::P][:V_bin]
                vmu1 = vmom[1::P][:V_bin] if P >= 2 else None
                vmu2 = vmom[2::P][:V_bin] if P >= 3 else None
                c_v = reconstruct_distribution(sf, vmu0, vmu1, vmu2,
                                               re.vac_bins, V)
                c_v[:v_d] = yj[iv:iv + v_d]
            else:
                # All discrete vacancies: only v_discrete entries stored
                c_v = np.zeros(V)
                c_v[:v_d] = yj[iv:iv + v_d]
        else:
            c_n = yj[:I]
            c_v = yj[I:I + V]
 
        # Fraction of content in last 10% of sizes
        tail_n = int(max(I * 0.1, 1))
        tail_m = int(max(V * 0.1, 1))
        total_sia = np.dot(ns, np.maximum(c_n, 0.0))
        total_vac = np.dot(ms, np.maximum(c_v, 0.0))
        tail_sia  = np.dot(ns[-tail_n:], np.maximum(c_n[-tail_n:], 0.0))
        tail_vac  = np.dot(ms[-tail_m:], np.maximum(c_v[-tail_m:], 0.0))
 
        # Guard: if total content is at floor level (all concentrations ≈ C_floor),
        # the tail fraction is a floor artifact, not real boundary pileup.
        # Require total content to be above I * C_floor * I (roughly I^2 * C_floor)
        # to distinguish real content from uniform-floor noise.
        C_floor = float(self.input_data.reactions.get('C_floor', 1e-15))
        sia_threshold = I * I * C_floor * 10.0
        vac_threshold = V * V * C_floor * 10.0
        frac_I = tail_sia / max(total_sia, 1e-300) if total_sia > sia_threshold else 0.0
        frac_V = tail_vac / max(total_vac, 1e-300) if total_vac > vac_threshold else 0.0
        return frac_I, frac_V
 
    def _boundary_fraction(self, results):
        """Check tail fractions at the last time step (backward compat)."""
        y = results.get('y')
        if y is None:
            return 0.0, 0.0
        return self._boundary_fraction_at(results, -1)
 
    def _find_first_exceedance(self, results, threshold):
        """
        Scan all output time points and return the index of the first one
        where either the SIA or vacancy tail fraction exceeds *threshold*.
 
        Returns None if no exceedance is found.
        """
        y = results.get('y')
        if y is None:
            return None
        n_t = y.shape[1]
        for j in range(n_t):
            frac_I, frac_V = self._boundary_fraction_at(results, j)
            if frac_I > threshold or frac_V > threshold:
                return j
        return None
 
    @staticmethod
    def _sia100_column(results, idx):
        """The ⟨100⟩ block AS SOLVED at time index *idx*, or None if absent.
 
        Reads ``y_sia100_raw`` (the raw appended block) and never ``y_sia100``,
        which is a per-size closure reconstruction with the C_floor IC already
        subtracted — fine for plots, wrong as an initial condition.
        """
        blk = results.get('y_sia100_raw')
        if blk is None:
            return None
        blk = np.asarray(blk, dtype=float)
        if blk.ndim != 2 or blk.size == 0:
            return None
        return blk[:, idx]
 
    def _resume_state(self, results, idx=-1):
        """Full C++ state vector at time index *idx*, ready to resume from.
 
        ``results['y']`` has the appended ⟨100⟩ block split off by cpp_bridge,
        so resuming from it alone hands the solver a SHORT y0.  The C++ reads
        missing ``y0_k`` keys as 1e-100 (parameters.h), so this did not fail —
        it silently erased the entire ⟨100⟩ population at every segment
        boundary, which the population then had to regrow from nothing.  That
        produced the once-per-segment sawtooth in the ⟨100⟩ number density and
        a violent restart transient for CVODE.  Re-attach the block here.
        """
        y = np.asarray(results['y'], dtype=float)[:, idx]
        blk = self._sia100_column(results, idx)
        return y if blk is None else np.concatenate([y, blk])
 
    def _expand_state(self, y_old, I_new, V_new, y100_old=None):
        """
        Map ODE state from the current domain to a larger (I_new, V_new) domain.
 
        1. Reconstruct per-size distributions from old state.
        2. Pad with C_floor for new sizes beyond old I, V.
        3. Resize domain (rebuilds bins and rate equations).
        4. Project padded distributions onto new bins.
 
        Returns y0_new for the enlarged domain.
 
        Parameters
        ----------
        y100_old : ndarray or None
            The appended ⟨100⟩ block AS SOLVED (``results['y_sia100_raw']``),
            i.e. per-size in discrete modes and discrete-prefix + bin moments in
            bin-moment modes.  When given it is remapped by the SAME
            reconstruct → pad → project path as the ½⟨111⟩ block (the two share
            a layout) and appended to the returned state, so the C++ receives a
            y0 of the full width.  Omitting it returns the ⟨100⟩-free state,
            which the solver zero-fills to 1e-100 — that erases the population
            at the doubling, so pass it whenever loop conversion is on.
        """
        re_old = self.rate_equations
        is_bin = hasattr(re_old, 'bins')
        I_old  = self.input_data.I
        V_old  = self.input_data.V
        C_floor = float(self.input_data.reactions.get('C_floor', 1e-15))
 
        yj = np.maximum(y_old, 0.0)
 
        if is_bin:
            from .bin_moment_rates import (reconstruct_distribution,
                                           moments_from_distribution)
            K_i_old = re_old.K_i
            K_v_old = re_old.K_v
            iv_old  = re_old.i_VAC
            P_old   = getattr(re_old, 'n_mom', 2)
            sf_old  = getattr(re_old, 'shape_function', 'linear')
            i_d_old = getattr(re_old, 'i_discrete', 0)
            v_d_old = getattr(re_old, 'v_discrete', 0)
 
            # Reconstruct SIA per-size from moments
            sia_mom = yj[i_d_old:i_d_old + P_old * K_i_old]
            sia_mu0 = sia_mom[0::P_old][:K_i_old]
            sia_mu1 = sia_mom[1::P_old][:K_i_old] if P_old >= 2 else None
            sia_mu2 = sia_mom[2::P_old][:K_i_old] if P_old >= 3 else None
            c_n_old = reconstruct_distribution(sf_old, sia_mu0, sia_mu1,
                                               sia_mu2, re_old.bins, I_old)
            c_n_old[:i_d_old] = yj[:i_d_old]
 
            # ⟨100⟩ block — layout matches the ½⟨111⟩ block, but it must NOT go
            # through the reconstruct → project path used above.  That path
            # round-trips the moments through `reconstruct_distribution`, whose
            # hat closure clips its negative lobe to zero; the clip is what
            # breaks μ₀ conservation.  For the ½⟨111⟩ population the effect is
            # mild (broad distribution, n̄ near bin centres).  The ⟨100⟩
            # population is the opposite — sessile, one-way, sharply peaked and
            # advecting in size space — so n̄ sits far from the bin centre and
            # the clip is catastrophic: measured ×176.9 inflation of μ₀ on the
            # real state at 1.9 dpa of the I=40000 run (2026-08-15), which is
            # exactly the ×177 jump in N_100 seen across that run's doubling.
            # Remap its moments conservatively instead — no per-size detour.
            mu0_100 = mu1_100 = None
            y100j = None
            if y100_old is not None:
                y100j = np.maximum(np.asarray(y100_old, dtype=float), 0.0)
                m100 = y100j[i_d_old:i_d_old + P_old * K_i_old]
                mu0_100 = m100[0::P_old][:K_i_old].astype(float)
                mu1_100 = (m100[1::P_old][:K_i_old].astype(float)
                           if P_old >= 2 else None)
 
            # Reconstruct vacancy per-size from moments
            vac_mom = yj[iv_old + v_d_old:iv_old + v_d_old + P_old * K_v_old]
            vmu0 = vac_mom[0::P_old][:K_v_old]
            vmu1 = vac_mom[1::P_old][:K_v_old] if P_old >= 2 else None
            vmu2 = vac_mom[2::P_old][:K_v_old] if P_old >= 3 else None
            c_v_old = reconstruct_distribution(sf_old, vmu0, vmu1, vmu2,
                                               re_old.vac_bins, V_old)
            c_v_old[:v_d_old] = yj[iv_old:iv_old + v_d_old]
 
            # Extract He state
            if re_old.he_mode == 'case2':
                Q_tot = yj[re_old.i_Qtot]
                Q_k_old = None
                q_n_old = None
            else:
                Q_k_old = yj[re_old.i_Q:re_old.i_Q + K_v_old]
                Q_tot = np.sum(Q_k_old)
                # Case 1 tracks per-bin mean He loading.  Reconstruct the
                # per-size He distribution q_n = ell_bin_k * c_v_n from the
                # OLD bins (ell_bin_k = Q_k / mu0_k), so the per-bin mean He
                # loading is preserved rather than smeared across the new
                # domain proportionally to mu0.  Same reconstruct->project
                # pattern as the SIA/VAC moment blocks above.
                q_n_old = np.zeros(V_old)
                for k, (mlo, mhi) in enumerate(re_old.vac_bins):
                    mu0_k = vmu0[k]
                    if mu0_k <= 0.0:
                        continue
                    ell_bin = Q_k_old[k] / mu0_k
                    ns = np.arange(mlo, mhi)
                    idx = ns - 1
                    valid = (idx >= 0) & (idx < V_old)
                    q_n_old[idx[valid]] = ell_bin * c_v_old[idx[valid]]
 
            c_h = None if re_old.qss_He else yj[re_old.i_He]
 
            # Extend distributions to new domain (pad with C_floor)
            c_n_new = np.full(I_new, C_floor)
            c_n_new[:I_old] = c_n_old
            c_v_new = np.full(V_new, C_floor)
            c_v_new[:V_old] = c_v_old
            # Per-size He distribution padded with 0 (new sizes hold no He)
            if q_n_old is not None:
                q_n_new = np.zeros(V_new)
                q_n_new[:V_old] = q_n_old
            else:
                q_n_new = None
 
            # Resize domain — rebuilds bins and rate equations
            self._resize_domain(I_new, V_new)
            re_new = self.rate_equations
            P_new = re_new.n_mom
 
            # Project onto new bins
            sia_mu0_new, sia_mu1_new, sia_mu2_new = moments_from_distribution(
                c_n_new, re_new.bins, n_mom=P_new)
            vac_mu0_new, vac_mu1_new, vac_mu2_new = moments_from_distribution(
                c_v_new, re_new.vac_bins, n_mom=P_new)
 
            # Build new state vector
            y0 = np.full(re_new.N_eq, C_floor)
            # Discrete SIA
            i_d_new = re_new.i_discrete
            y0[:i_d_new] = c_n_new[:i_d_new]
            # Binned SIA moments
            for k in range(re_new.K_i):
                y0[i_d_new + P_new * k] = sia_mu0_new[k]
                if P_new >= 2 and sia_mu1_new is not None:
                    y0[i_d_new + P_new * k + 1] = sia_mu1_new[k]
                if P_new >= 3 and sia_mu2_new is not None:
                    y0[i_d_new + P_new * k + 2] = sia_mu2_new[k]
            # Discrete VAC
            iv_new = re_new.i_VAC
            v_d_new = re_new.v_discrete
            y0[iv_new:iv_new + v_d_new] = c_v_new[:v_d_new]
            # Binned VAC moments
            vac_start_new = iv_new + v_d_new
            for k in range(re_new.K_v):
                y0[vac_start_new + P_new * k] = vac_mu0_new[k]
                if P_new >= 2 and vac_mu1_new is not None:
                    y0[vac_start_new + P_new * k + 1] = vac_mu1_new[k]
                if P_new >= 3 and vac_mu2_new is not None:
                    y0[vac_start_new + P_new * k + 2] = vac_mu2_new[k]
 
            # He state
            if re_new.he_mode == 'case2':
                y0[re_new.i_Qtot] = Q_tot
            else:
                # Case 1: re-project the reconstructed per-size He
                # distribution q_n onto the NEW vacancy bins.  Q_k for a
                # bin is exactly Σ_{n∈bin} q_n, i.e. the zeroth moment of
                # q_n over that bin — this preserves the per-bin mean He
                # loading instead of smearing Q_tot across bins by mu0.
                if q_n_new is not None:
                    Q_k_new, _, _ = moments_from_distribution(
                        q_n_new, re_new.vac_bins, n_mom=1)
                    for k in range(re_new.K_v):
                        y0[re_new.i_Q + k] = Q_k_new[k]
            if not re_new.qss_He and c_h is not None:
                y0[re_new.i_He] = c_h
 
            # Carry the 5 conservation accounting integrals (J_SIA_fixed,
            # J_SIA_mutual, J_VAC_fixed, J_VAC_mutual, J_He_sink — last 5
            # state entries) across the doubling; resetting them to C_floor
            # corrupts delta_FP / delta_He for the remainder of the run.
            y0[-5:] = yj[-5:]
 
            if mu0_100 is None:
                return y0
            # Conservative moment remap (see the note where mu0_100 is read).
            mu0_n, mu1_n = _remap_bin_moments(mu0_100, mu1_100,
                                              re_old.bins, re_new.bins)
            blk = np.full(i_d_new + P_new * re_new.K_i, C_floor)
            # i_discrete does not change on a doubling, so the discrete prefix
            # carries over one-for-one; guard the general case anyway.
            n_pre = min(i_d_new, i_d_old, y100j.size)
            blk[:n_pre] = y100j[:n_pre]
            for k in range(re_new.K_i):
                blk[i_d_new + P_new * k] = mu0_n[k]
                if P_new >= 2:
                    blk[i_d_new + P_new * k + 1] = mu1_n[k]
                if P_new >= 3:
                    # No independent μ₂ transport in the remap: place it
                    # consistently with the bin's own mean so the lognormal
                    # closure sees a self-consistent triple rather than a
                    # μ₂ inherited from a different binning.
                    nb = (mu1_n[k] / mu0_n[k]) if mu0_n[k] > 0 else 0.0
                    blk[i_d_new + P_new * k + 2] = mu0_n[k] * nb * nb
            return np.concatenate([y0, blk])
 
        else:
            # Full per-size mode — pad state vector directly
            c_n_old = yj[:I_old]
            c_v_old = yj[I_old:I_old + V_old]
 
            if re_old.he_mode == 'case2':
                Q_tot = yj[re_old.i_Qtot]
                Q_m_old = None
            else:
                Q_m_old = yj[re_old.i_Q:re_old.i_Q + V_old]
                Q_tot = np.sum(Q_m_old)
            c_h = None if re_old.qss_He else yj[re_old.i_He]
 
            # Resize domain
            self._resize_domain(I_new, V_new)
            re_new = self.rate_equations
 
            y0 = np.full(re_new.N_eq, C_floor)
            y0[:I_old] = c_n_old
            y0[I_new:I_new + V_old] = c_v_old
 
            if re_new.he_mode == 'case2':
                y0[re_new.i_Qtot] = Q_tot
            else:
                y0[re_new.i_Q:re_new.i_Q + V_old] = Q_m_old
            if not re_new.qss_He and c_h is not None:
                y0[re_new.i_He] = c_h
 
            # Carry the 5 conservation accounting integrals across the
            # doubling (see bin-moment branch above).
            y0[-5:] = yj[-5:]
 
            if y100_old is None:
                return y0
            # Discrete mode: the ⟨100⟩ block is already per-size, so widening it
            # is a pad, exactly as for the ½⟨111⟩ block.
            y100j = np.maximum(np.asarray(y100_old, dtype=float), 0.0)
            blk = np.full(I_new, C_floor)
            blk[:min(I_old, y100j.size)] = y100j[:min(I_old, y100j.size)]
            return np.concatenate([y0, blk])
 
    # ── Time-series keys that should be concatenated when merging segments ────
 
    _TS_KEYS = {'t', 'dose', 'rho_net', 'C_SIA_tot', 'C_VAC_tot', 'C_He_tot',
                'C_He_free', 'mean_n_i', 'mean_n_v', 'N_loops', 'N_voids',
                'swelling', 'C_i1', 'C_v1', 'delta_FP', 'delta_He',
                'J_SIA_fixed', 'J_SIA_mutual', 'J_VAC_fixed', 'J_VAC_mutual',
                'J_He_sink',
                # loop-conversion series (⟨100⟩ block); omitting these left
                # them holding only the LAST segment, so their length no
                # longer matched len(t) and every loop-character plot broke.
                'N_loops_111', 'N_loops_100', 'mean_n_111', 'mean_n_100',
                'f_111_loop',
                # the separated conservation arms — same defect
                'delta_FP_sia', 'delta_FP_vac'}
 
    # 2-D [size, time] arrays: sliced by column, merged with row zero-padding
    # so a domain doubling (which grows the size axis) still concatenates.
    _TS_KEYS_2D = {'y_sia100', 'y_sia100_raw'}
 
    def _slice_results(self, results, start, end):
        """Slice a results dict to time indices [start:end)."""
        sliced = {}
        n_t = len(results['t'])
        for key, val in results.items():
            if key in self._TS_KEYS and isinstance(val, np.ndarray):
                sliced[key] = val[start:end]
            elif (key == 'y' or key in self._TS_KEYS_2D) and \
                    isinstance(val, np.ndarray) and val.ndim == 2:
                sliced[key] = val[:, start:end]
            else:
                sliced[key] = val
        return sliced
 
    def _merge_results(self, accumulated, new_segment):
        """
        Merge two results dicts by concatenating time-series arrays.
 
        The first time point of *new_segment* is skipped if it duplicates
        the last time point of *accumulated* (continuation overlap).
        The raw state ``y`` is decomposed into SIA / VAC / He blocks,
        each zero-padded independently to the larger layout, then
        reassembled and concatenated so that plots cover the full dose
        range even after domain-doubling.
        """
        if accumulated is None:
            return new_segment
        if new_segment is None:
            return accumulated
 
        # Detect overlap: first point of new matches last of accumulated
        skip = 0
        if (len(accumulated['t']) > 0 and len(new_segment['t']) > 0 and
                abs(new_segment['t'][0] - accumulated['t'][-1]) /
                max(abs(accumulated['t'][-1]), 1e-30) < 1e-8):
            skip = 1
 
        merged = {}
        for key in new_segment:
            old_val = accumulated.get(key)
            new_val = new_segment[key]
            if key in self._TS_KEYS and isinstance(new_val, np.ndarray):
                if old_val is not None and isinstance(old_val, np.ndarray):
                    merged[key] = np.concatenate([old_val, new_val[skip:]])
                else:
                    merged[key] = new_val
            elif key == 'y':
                if old_val is not None and isinstance(old_val, np.ndarray):
                    merged[key] = self._merge_y_blocks(
                        accumulated, new_segment, skip)
                else:
                    merged[key] = new_val
            elif key in self._TS_KEYS_2D and isinstance(new_val, np.ndarray) \
                    and new_val.ndim == 2:
                if old_val is not None and isinstance(old_val, np.ndarray) \
                        and old_val.ndim == 2:
                    merged[key] = self._concat_2d_padded(
                        old_val, new_val[:, skip:])
                else:
                    merged[key] = new_val
            else:
                merged[key] = new_val
        return merged
 
    @staticmethod
    def _concat_2d_padded(old, new):
        """Concatenate two [size, time] arrays along time, zero-padding the
        shorter size axis so a domain doubling does not break the merge."""
        n_old, n_new = old.shape[0], new.shape[0]
        if n_old < n_new:
            old = np.vstack([old, np.zeros((n_new - n_old, old.shape[1]))])
        elif n_new < n_old:
            new = np.vstack([new, np.zeros((n_old - n_new, new.shape[1]))])
        return np.concatenate([old, new], axis=1)
 
    @staticmethod
    def _merge_y_blocks(accumulated, new_segment, skip):
        """
        Remap old y to the new layout and concatenate.
 
        Uses ``_y_i_VAC`` and ``_y_i_He`` stored in each results dict to
        decompose y into SIA / VAC / He blocks, zero-pad each to the
        larger size, then reassemble into a single array with the new
        layout.
        """
        old_y  = accumulated['y']
        new_y  = new_segment['y']
 
        iv_old = accumulated.get('_y_i_VAC', old_y.shape[0])
        ih_old = accumulated.get('_y_i_He',  old_y.shape[0])
        iv_new = new_segment.get('_y_i_VAC', new_y.shape[0])
        ih_new = new_segment.get('_y_i_He',  new_y.shape[0])
 
        # The SIA/VAC/He layout indices are required to remap old y into
        # the new layout.  Silently returning only the latest segment here
        # discards accumulated history and yields a `y` whose column count
        # no longer matches len(t) — fail loudly instead.
        if '_y_i_VAC' not in accumulated or '_y_i_VAC' not in new_segment:
            missing = ('accumulated' if '_y_i_VAC' not in accumulated
                       else 'new_segment')
            raise ValueError(
                "_merge_y_blocks: '_y_i_VAC' missing from "
                f"{missing} results dict — cannot remap y blocks without "
                "the SIA/VAC/He layout indices. The merged y would have a "
                "column count inconsistent with len(t).")
 
        # Decompose into blocks
        sia_old = old_y[:iv_old, :]
        vac_old = old_y[iv_old:ih_old, :]
        he_old  = old_y[ih_old:, :]
 
        sia_new = new_y[:iv_new, :]
        vac_new = new_y[iv_new:ih_new, :]
        he_new  = new_y[ih_new:, :]
 
        # Zero-pad each block of old to match new sizes
        def _pad_rows(arr, target_rows):
            if arr.shape[0] >= target_rows:
                return arr[:target_rows, :]
            pad = np.zeros((target_rows - arr.shape[0], arr.shape[1]))
            return np.vstack([arr, pad])
 
        sia_old = _pad_rows(sia_old, iv_new)
        vac_old = _pad_rows(vac_old, ih_new - iv_new)
        he_n    = he_new.shape[0]
        he_old  = _pad_rows(he_old, he_n)
 
        # Reassemble old y in new layout
        old_remapped = np.vstack([sia_old, vac_old, he_old])
        return np.concatenate([old_remapped, new_y[:, skip:]], axis=1)
 
    def run_adaptive(self, solver_config=None, save_output=True,
                     progress_callback=None, timeout_s=None,
                     boundary_threshold=0.01, max_doublings=6,
                     points_per_segment=10):
        """
        Truly adaptive domain doubling via short time-segments.
 
        The total time span is divided into short segments of
        *points_per_segment* output points each.  After every segment
        the tail fraction is checked.  If it exceeds *boundary_threshold*:
 
        1. The first exceedance point within the segment is located.
        2. Results up to that point are kept.
        3. I and/or V are doubled; the state vector is mapped to the
           larger domain via ``_expand_state``.
        4. Integration resumes from the exceedance time.
 
        This allows starting with a small domain (e.g. I = V = 1000) and
        growing it on demand as the distribution spreads, without ever
        restarting from t = 0.
 
        Parameters
        ----------
        boundary_threshold : float
            Adapt if tail fraction > this (default 0.01 = 1 %).
        max_doublings : int
            Maximum number of doublings (default 6; 1000 → 64 000).
        points_per_segment : int
            Output points per solver call (default 10).  Smaller values
            detect exceedance earlier but add subprocess overhead.
 
        Returns
        -------
        results : dict
        """
        if solver_config is None:
            solver_config = self._default_solver_config()
 
        t_begin, t_end = solver_config['t_span']
        n_points_total = solver_config.get('n_points', 100)
        log_time       = solver_config.get('log_time', True)
 
        # Precompute checkpoint times (segment boundaries)
        n_segments = max(1, n_points_total // points_per_segment)
        if log_time and t_begin > 0:
            checkpoints = np.geomspace(t_begin, t_end, n_segments + 1)
        else:
            checkpoints = np.linspace(t_begin, t_end, n_segments + 1)
 
        # REQUIRED TIMES (added 2026-09-06).  The output grid is logarithmic in
        # TIME, so the doses a study reports are not generally on it.  Measured:
        # a 40 dpa row with n_points=37 and points_per_segment=10 saves 28
        # points, of which only four lie above 1 dpa -- 0.954, 3.315, 11.51, 40.
        # A request to score at 15.72 dpa then silently took 11.51 dpa, 26.8%
        # low, because the state at 15.72 was never written.  Nothing downstream
        # could repair that: it is missing data, not a lookup error.
        #
        # Making each required time a SEGMENT BOUNDARY puts it in the output
        # exactly, because every segment is integrated to its own endpoint.
        # That is why this belongs here and not in a post-hoc resample: the
        # solver has to stop there.  Everything else is unchanged -- the log
        # spacing, points_per_segment, and the segment machinery all still
        # apply, with a few extra boundaries in the range that matters.
        req = solver_config.get('required_times') or ()
        req = [float(x) for x in req
               if t_begin < float(x) < t_end * (1 - 1e-12)]
        if req:
            checkpoints = np.unique(np.concatenate([checkpoints, req]))
            # n_segments is now len-1, and max_iters below is derived from it;
            # leaving the old value would cap the loop before the last segment.
            n_segments = len(checkpoints) - 1
 
        # ROW-LEVEL DEADLINE, not a per-segment one.  timeout_s used to be handed
        # to EVERY segment unchanged, so a run could legally take
        # n_segments x timeout_s: measured 2026-08-06, rows of 16108-20503 s
        # completed under a 12000 s setting, and the "time limit" the campaign
        # was designed around never bit.  Budget it across the whole run instead
        # and let the last segment be cut gracefully -- the solver flushes what
        # it has, so the row still contributes every dose-ladder rung it reached
        # (plan S12(s)).
        import time as _time
        _deadline = (_time.monotonic() + timeout_s) if timeout_s else None
 
        accumulated  = None
        self._accumulated_results = None   # expose for graceful interrupt
        current_t    = t_begin
        y0_override  = None
        n_doublings_I = 0     # independent doubling counters
        n_doublings_V = 0
        seg_count    = 0
        interrupted  = False
        _deferred_exc = None   # a crash is re-raised AFTER the save block
        # Allow extra iterations for expansion restarts within segments
        max_iters    = n_segments + 2 * max_doublings + 1
 
        try:
          while current_t < t_end * (1 - 1e-10) and seg_count < max_iters:
            seg_count += 1
            I_cur = self.input_data.I
            V_cur = self.input_data.V
 
            # Next checkpoint after current_t
            future = checkpoints[checkpoints > current_t * (1 + 1e-10)]
            seg_t_end = float(future[0]) if len(future) > 0 else t_end
 
            print(f"\n[Segment {seg_count}] I={I_cur}  V={V_cur}"
                  f"  t=[{current_t:.2e}, {seg_t_end:.2e}]")
 
            seg_config = dict(solver_config)
            seg_config['t_span']   = (current_t, seg_t_end)
            seg_config['n_points'] = points_per_segment
 
            # Whatever is LEFT of the row's budget, not the full budget again.
            if _deadline is not None:
                _left = _deadline - _time.monotonic()
                if _left <= 0:
                    print(f"  Row budget of {timeout_s:.0f}s exhausted after "
                          f"{seg_count-1} segment(s) — stopping with what is "
                          f"integrated so far.")
                    interrupted = True
                    break
            else:
                _left = None
 
            results = self.run(
                solver_config=seg_config,
                save_output=False,
                progress_callback=progress_callback,
                timeout_s=_left,
                y0_override=y0_override,
            )
 
            if results is None:
                print("  Solver failed — returning accumulated results.")
                if accumulated is not None and save_output:
                    self._diag_text = self.reaction_rates.format_diagnostic(
                        mean_n_i=accumulated['mean_n_i'][-1]
                        if 'mean_n_i' in accumulated else None)
                    self._save_output(accumulated, solver_config)
                return accumulated
 
            # If this segment came back partial (Ctrl+C / timeout / signal),
            # honor the interrupt: merge what we have and stop the adaptive
            # loop instead of launching another segment.
            seg_partial = results.get('metadata', {}).get(
                'solver_stats', {}).get('partial', False)
            if seg_partial:
                print("  Segment ended early (interrupt or timeout) — "
                      "stopping adaptive run.")
                accumulated = self._merge_results(accumulated, results)
                self._accumulated_results = accumulated
                interrupted = True
                break
 
            # Operator-split advance of the dynamic network density ρ_net, once
            # per segment regardless of whether the domain doubles below (the
            # doubling branch does not reach the else-clause).  No-op unless
            # LOOP_NETWORK_LOSS is on.  Freezes the segment-end monomers and
            # rebuilds the rates so the next segment picks up the new ρ_net.
            self._loop_network_update(
                results,
                float(results['t'][-1] - results['t'][0]))
 
            # Check tail at the last point of this segment
            frac_I, frac_V = self._boundary_fraction_at(results, -1)
            can_double_I = frac_I > boundary_threshold and n_doublings_I < max_doublings
            can_double_V = frac_V > boundary_threshold and n_doublings_V < max_doublings
 
            if can_double_I or can_double_V:
                # Find the first exceedance point within this segment
                exceed_idx = self._find_first_exceedance(
                    results, boundary_threshold)
                if exceed_idx is None:
                    exceed_idx = results['y'].shape[1] - 1
 
                frac_I_ex, frac_V_ex = self._boundary_fraction_at(
                    results, exceed_idx)
                t_exceed = results['t'][exceed_idx]
 
                # Keep results up to and including exceedance
                partial = self._slice_results(results, 0, exceed_idx + 1)
                accumulated = self._merge_results(accumulated, partial)
                self._accumulated_results = accumulated
 
                # Double each dimension independently if it exceeds AND
                # still has budget remaining
                do_I = frac_I_ex > boundary_threshold and n_doublings_I < max_doublings
                do_V = frac_V_ex > boundary_threshold and n_doublings_V < max_doublings
                I_new = I_cur * 2 if do_I else I_cur
                V_new = V_cur * 2 if do_V else V_cur
                which = []
                if do_I:
                    which.append(f'I: {I_cur}→{I_new} [{n_doublings_I+1}/{max_doublings}]')
                    n_doublings_I += 1
                if do_V:
                    which.append(f'V: {V_cur}→{V_new} [{n_doublings_V+1}/{max_doublings}]')
                    n_doublings_V += 1
                print(f"  Exceedance at t={t_exceed:.3e}: "
                      f"SIA={frac_I_ex:.3f}  VAC={frac_V_ex:.3f}")
                print(f"  Doubling: {', '.join(which)}")
 
                # Map state to enlarged domain.  The ⟨100⟩ block rides along —
                # see _resume_state below for why it must.
                y_at   = results['y'][:, exceed_idx]
                y100_at = self._sia100_column(results, exceed_idx)
                y0_override = self._expand_state(y_at, I_new, V_new, y100_at)
                current_t   = t_exceed
            else:
                # Segment OK (or both dimensions maxed out) — accumulate
                exceeded_but_maxed = []
                if frac_I > boundary_threshold and n_doublings_I >= max_doublings:
                    exceeded_but_maxed.append(f'SIA={frac_I:.3f}')
                if frac_V > boundary_threshold and n_doublings_V >= max_doublings:
                    exceeded_but_maxed.append(f'VAC={frac_V:.3f}')
                if exceeded_but_maxed:
                    print(f"  Tail: {', '.join(exceeded_but_maxed)} "
                          f"(max doublings reached)")
                else:
                    print(f"  Tail OK: SIA={frac_I:.3f}  VAC={frac_V:.3f}")
 
                accumulated = self._merge_results(accumulated, results)
                self._accumulated_results = accumulated
                current_t   = results['t'][-1]
                y0_override = self._resume_state(results, -1)
 
        except KeyboardInterrupt:
            interrupted = True
            print("\n\n*** KeyboardInterrupt — stopping gracefully. ***")
            # Merge the last completed segment if it wasn't merged yet
            if results is not None and results is not accumulated:
                try:
                    accumulated = self._merge_results(accumulated, results)
                    self._accumulated_results = accumulated
                except Exception:
                    pass   # keep whatever was accumulated before
        except BaseException as _exc:
            # Same reasoning as the orchestrator: a crash in segment k must not
            # discard segments 1..k-1.  A long adaptive run is exactly where
            # this costs most -- the accumulated trajectory can represent many
            # hours.  Merge whatever the failing segment returned, fall through
            # to the normal save block below, and re-raise afterwards so the
            # error is still reported.
            interrupted = True
            _deferred_exc = _exc
            print(f"\n\n*** {type(_exc).__name__} in segment {seg_count} — "
                  f"saving the {seg_count - 1} completed segment(s) before "
                  f"re-raising. ***")
            if results is not None and results is not accumulated:
                try:
                    accumulated = self._merge_results(accumulated, results)
                    self._accumulated_results = accumulated
                except Exception:
                    pass
 
        # ── Finished (or interrupted) — save and return ───────────────────
        if accumulated is not None:
            I_final = self.input_data.I
            V_final = self.input_data.V
            tot = n_doublings_I + n_doublings_V
            status = "INTERRUPTED" if interrupted else "complete"
            n_pts  = len(accumulated['t'])
            print(f"\nAdaptive run {status}: {seg_count} segments, "
                  f"{tot} doublings (I×{n_doublings_I}, V×{n_doublings_V}), "
                  f"final domain I={I_final} V={V_final}, "
                  f"{n_pts} time points saved")
            if save_output:
                self._diag_text = self.reaction_rates.format_diagnostic(
                    mean_n_i=accumulated['mean_n_i'][-1]
                    if 'mean_n_i' in accumulated else None)
                self._save_output(accumulated, solver_config)
        elif interrupted:
            print("\nNo completed segments to save.")
 
        # Re-raise a crash only now, so the traceback still reaches the caller
        # but the completed segments are already on disk.  A KeyboardInterrupt
        # deliberately does NOT come through here: stopping early is a normal
        # outcome for this code and returning the partial trajectory is the
        # documented behaviour.
        if _deferred_exc is not None:
            raise _deferred_exc
 
        return accumulated
 
    # ── Main entry point ──────────────────────────────────────────────────────
 
    def run(self, solver_config=None, save_output=True, progress_callback=None,
            timeout_s=None, y0_override=None):
        """
        Run the simulation using the configured solver mode and physics option.
 
        Parameters
        ----------
        solver_config     : dict, optional
            Keys: t_span, n_points, rtol, atol, log_time,
                  solver_method (dict with linsol, window_*, etc.)
        save_output       : bool
            Write timestamped output/ directory.
        progress_callback : callable or None
            Called once per output time step during the C++ solver run with a
            dict of concentrations and rate-breakdown values (atom fraction).
            Enables verbose mode in the C++ solver automatically.
        timeout_s         : float or None
            Maximum wall-clock seconds for the C++ solver.
        y0_override       : ndarray or None
            Custom initial conditions (adaptive continuation).
 
        Returns
        -------
        results : dict  (see post_process.calculate_derived_quantities)
        """
        if solver_config is None:
            solver_config = self._default_solver_config()
 
        # Explicit kwarg wins; otherwise fall back to the config-level default
        # (added by _default_solver_config: 1 hour wall-clock cap).
        if timeout_s is None:
            timeout_s = solver_config.get('timeout_s')
 
        sm = self.input_data.solver_mode
        print(f"\nLaunching solver_mode='{sm}' …")
 
        _t0 = _time.perf_counter()
        results = None
        # Pre-clear any stash from a previous run so we don't accidentally
        # save stale data from an earlier interrupted run.
        self._partial_results = None
        try:
            if sm in ('full_system', 'active_window'):
                results = self._run_cpp(solver_config, progress_callback,
                                        timeout_s=timeout_s,
                                        y0_override=y0_override)
            else:
                raise ValueError(f"Unknown solver_mode='{sm}'. "
                                 "Use full_system or active_window.")
        except KeyboardInterrupt:
            # Defense-in-depth: cpp_bridge already catches Ctrl+C and returns
            # whatever it could parse, but a *second* Ctrl+C during its own
            # graceful-shutdown wait can propagate up to here.  Recover from
            # the stash that cpp_bridge sets as soon as it has any data, so
            # the partial output still gets saved.
            print("\n*** KeyboardInterrupt at orchestrator — "
                  "recovering partial output. ***")
            results = getattr(self, '_partial_results', None)
        except BaseException as _exc:
            # A CRASH MUST NOT COST THE TRAJECTORY.  Ctrl+C was already handled
            # above, but any other failure -- solver segfault surfacing as an
            # exception, MemoryError, a parse error on a truncated stdout --
            # previously propagated with nothing written, discarding every hour
            # already spent.  Recover cpp_bridge's stash, flush it to a
            # timestamped directory exactly as a normal run would, and only then
            # re-raise so the traceback still reaches the user.
            _partial = getattr(self, '_partial_results', None)
            print(f"\n*** {type(_exc).__name__} during solve — "
                  f"{'flushing partial output' if _partial is not None else 'no partial output to save'}"
                  f" before re-raising. ***")
            if _partial is not None:
                self._accumulated_results = _partial
                self._wall_clock_s = _time.perf_counter() - _t0
                try:
                    self._diag_text = self.reaction_rates.format_diagnostic(
                        mean_n_i=_partial['mean_n_i'][-1]
                        if 'mean_n_i' in _partial else None)
                except Exception:
                    self._diag_text = ''
                if save_output:
                    try:
                        self._save_output(_partial, solver_config)
                    except Exception as _se:      # never mask the real error
                        print(f"    (save failed: {type(_se).__name__}: {_se})")
            raise
        self._wall_clock_s = _time.perf_counter() - _t0
 
        # Mirror the cpp_bridge stash to the run_adaptive-style attribute so
        # the notebook fallback `if results is None: results = sim._accumulated_results`
        # works for the non-adaptive path too.
        if results is not None:
            self._accumulated_results = results
            # Store diagnostic text for file output (no inline printing)
            try:
                self._diag_text = self.reaction_rates.format_diagnostic(
                    mean_n_i=results['mean_n_i'][-1] if 'mean_n_i' in results else None
                )
            except Exception:
                self._diag_text = ''
            if save_output:
                self._save_output(results, solver_config)
 
        return results
 
    # ── C++ solver dispatch ───────────────────────────────────────────────────
 
    def _run_cpp(self, solver_config, progress_callback=None, timeout_s=None,
                 y0_override=None):
        """Invoke the C++ solver via cpp_bridge."""
        from . import cpp_bridge
        results = cpp_bridge.run_cpp_solver(
            self, solver_config, base_dir=BASE_DIR,
            progress_callback=progress_callback,
            timeout_s=timeout_s,
            y0_override=y0_override,
        )
        return results
 
    # ── Default solver config ─────────────────────────────────────────────────
 
    def _default_solver_config(self):
        re = self.input_data.reactions
        sm = self.input_data.solver_mode
        po = self.input_data.physics_option
 
        # Linear solver and window parameters from reactions sheet
        linsol = str(re.get('linsol', 'dense')).lower()
        w_w    = int(float(re.get('window_width', 100)))
        C_exp  = float(re.get('window_C_exp',  1e-18))
 
        # Map solver mode to window_mode integer
        window_mode_map = {
            'full_system':   0,
            'active_window': 4,
        }
        win_mode = window_mode_map.get(sm, 0)
 
        # For bin_moment, force gmres: dense/band scale poorly at large N_eq
        # and klu's sparse direct path is only wired up for full_CD.
        if 'bin_moment' in po:
            linsol = 'gmres'
 
        t_begin = float(re.get('t_begin', 1e-8))
        t_end   = float(re.get('t_end',   1e7))
 
        return {
            't_span':    (t_begin, t_end),
            'n_points':  int(float(re.get('n_points', 200))),
            'log_time':  bool(int(float(re.get('log_time', 1)))),
            'rtol':      float(re.get('rtol', 1e-8)),
            'atol':      float(re.get('atol', 1e-20)),
            # Wall-clock cap for the C++ solver; partial output is saved on
            # timeout.  Override per-call via run(timeout_s=…) or by editing
            # this entry in solver_config.
            'timeout_s': float(re.get('timeout_s', 3600.0)),
            'solver_method': {
                'linsol':               linsol,
                'window_mode':          win_mode,
                'window_width':         w_w,
                'concentration_threshold': C_exp,
                'window_pad':              10,
                # CVODE max internal step.  Capping at t_end/10 prevents BDF
                # from speculating across many simulated decades in one step
                # (which leads to Newton failures at near-steady states).
                # 0 means no limit; override via solver_method['hmax'].
                'hmax':                 float(re.get('hmax', t_end / 10.0)),
            },
        }
 
    # ── Output writing ────────────────────────────────────────────────────────
 
    @staticmethod
    def _fmt_value(v):
        """Format a scalar for the provenance file."""
        if callable(v):
            return '<callable>'
        if isinstance(v, float):
            return f"{v:.6g}"
        if isinstance(v, (tuple, list)):
            return '[' + ', '.join(RadClusterSimulation._fmt_value(x) for x in v) + ']'
        return str(v)
 
    @classmethod
    def _write_dict_block(cls, f, title, d):
        f.write(f"### {title}\n\n")
        if not d:
            f.write("_(empty)_\n\n")
            return
        for k in sorted(d.keys(), key=str):
            f.write(f"- {k}: {cls._fmt_value(d[k])}\n")
        f.write("\n")
 
    def _machine_info(self):
        info = {
            'hostname':       platform.node(),
            'platform':       platform.platform(),
            'processor':      platform.processor() or platform.machine(),
            'python':         sys.version.split()[0],
            'cpu_count':      os.cpu_count(),
            'omp_num_threads': os.environ.get('OMP_NUM_THREADS', '<unset>'),
        }
        try:
            import psutil
            vm = psutil.virtual_memory()
            info['ram_total_GB'] = round(vm.total / 1024**3, 2)
            info['ram_available_GB'] = round(vm.available / 1024**3, 2)
            proc = psutil.Process(os.getpid())
            info['process_rss_GB'] = round(proc.memory_info().rss / 1024**3, 3)
        except Exception:
            info['ram_total_GB'] = '<psutil unavailable>'
        return info
 
    def _save_output(self, results, solver_config):
        """Write timestamped output directory.
 
        Resilient to partial/interrupted runs: each artifact (provenance,
        summary, diagnostics, .npy, plots) is written under its own try/except
        so a failure in one does not abort the others.
 
        The whole body runs under an `_InterruptDeferred` shield so that a
        Ctrl+C arriving mid-save does not abandon partially-written files —
        the current write completes, the rest of the artifacts are flushed,
        then the interrupt is re-raised for the caller.  A second Ctrl+C in
        the same block restores the default handler so the user can still
        force-abort if a save step hangs.
        """
        with _InterruptDeferred():
            return self._save_output_impl(results, solver_config)
 
    def _save_output_impl(self, results, solver_config):
        sm = self.input_data.solver_mode
        po = self.input_data.physics_option
        I  = self.input_data.I
        V  = self.input_data.V
        im = self.input_data.i_mobile
        vm = self.input_data.v_mobile
        ts = datetime.now().strftime('%Y%m%d_%H%M%S')
        domain   = f"I{I}V{V}"
        mobility = f"im{im}vm{vm}"
        # OPTIONAL CALLER TAG.  The label below is built from solver mode,
        # physics option, domain and mobility -- none of which distinguish the
        # rungs of a closure-refinement ladder, which differ only in
        # i_discrete/I_bin/v_discrete/V_bin.  Without a tag the four Table-1
        # rungs of the verification study produce four directory names that
        # differ only by timestamp.  Set `sim.run_tag = 'T1_B3'` before the run.
        tag = ''
        raw_tag = getattr(self, 'run_tag', None)
        if raw_tag:
            # Anything that could create a path level or confuse the label
            # parser is collapsed; the tag is a directory-name component.
            tag = '_' + re.sub(r'[^A-Za-z0-9._+-]+', '-', str(raw_tag)).strip('-_')[:48]
        # Tag interrupted runs in the directory name for easy identification.
        partial_tag = ''
        try:
            if results.get('metadata', {}).get('solver_stats', {}).get('partial'):
                partial_tag = '_PARTIAL'
        except Exception:
            pass
        # Tag sits directly after the timestamp: name order still equals time
        # order, and the rung is the first thing the eye reaches in a listing
        # of two dozen otherwise-identical run directories.
        label = f"{ts}{tag}_{sm}_{po}_{domain}_{mobility}{partial_tag}"
 
        out_dir  = BASE_DIR / 'output' / label
        plot_dir = out_dir / 'plots'
        plot_dir.mkdir(parents=True, exist_ok=True)
        print(f"Saving output to: {out_dir}")
        n_pts_have = len(results.get('t', [])) if isinstance(results, dict) else 0
        print(f"  {n_pts_have} time points to save")
 
        # ── Provenance (4-section layout) ─────────────────────────────────
        idata   = self.input_data
        derived = {k: v for k, v in idata.derived.items() if not callable(v)}
 
        # User selections — what the caller chose (post-override)
        user_selections = {
            'solver_mode':    sm,
            'physics_option': po,
            'I':              I,
            'V':              V,
            'i_mobile':       im,
            'v_mobile':       vm,
            'i_discrete':     idata.i_discrete,
            'v_discrete':     idata.v_discrete,
            'I_bin':          idata.I_bin,
            'V_bin':          idata.V_bin,
            'shape_function': idata.shape_function,
            'L_He_max':       idata.L_He_max,
            'alpha_He':       idata.alpha_He,
            'he_kinetics':    idata.reactions.get('he_kinetics', 'dynamic'),
            'C_floor':        idata.reactions.get('C_floor', 1e-15),
            'boundary_flux':  derived.get('boundary_flux'),
            'spectrum':       derived.get('spectrum'),
            'T':              derived.get('T'),
            'G':              derived.get('G'),
            'excel_file':     str(idata.excel_file.resolve()),
        }
 
        # Solver configuration — flatten nested solver_method
        sc_flat = {k: v for k, v in solver_config.items() if k != 'solver_method'}
        for k, v in solver_config.get('solver_method', {}).items():
            sc_flat[f"solver_method.{k}"] = v
 
        # Run statistics
        n_pts = results.get('metadata', {}).get('solver_stats', {}).get(
            'n_time_points', len(results.get('t', [])))
        run_stats = {
            'wall_clock_s':   getattr(self, '_wall_clock_s', None),
            'n_time_points':  n_pts,
            'timestamp':      ts,
        }
        run_stats.update(self._machine_info())
        omp_used = results.get('metadata', {}).get('omp_threads_used')
        if omp_used is not None:
            run_stats['omp_threads_used'] = omp_used
        ssf = results.get('metadata', {}).get('solver_stats_final')
        if isinstance(ssf, dict):
            for k in ('steps', 'nfe', 'nni', 'nli', 'nli_per_nni',
                      'npe', 'nps', 'ncfn', 'netf', 'nlsetup'):
                if k in ssf:
                    run_stats[f'solver.{k}'] = ssf[k]
        # Mark interrupted runs explicitly so partial output is identifiable
        # in the provenance without having to compare time-point counts.
        partial = results.get('metadata', {}).get('solver_stats', {}).get('partial')
        if partial is not None:
            run_stats['partial'] = bool(partial)
            run_stats['run_status'] = 'interrupted' if partial else 'completed'
 
        # ── 1. Provenance (always safe — pulls from self/input_data only) ────
        try:
            with open(out_dir / 'provenance.md', 'w', encoding='utf-8') as f:
                f.write(f"# RadCluster_2_1 run — {label}\n\n")
 
                f.write("## (1) Material Data\n\n")
                f.write("_All input tables from the Excel workbook with parameter "
                        "overrides applied._\n\n")
                self._write_dict_block(f, 'Production (Fission)',  idata.production_fission)
                self._write_dict_block(f, 'Production (Fusion)',   idata.production_fusion)
                self._write_dict_block(f, 'Energetics',            idata.energetics)
                self._write_dict_block(f, 'Diffusion',             idata.diffusion)
                self._write_dict_block(f, 'Dissociation',          idata.dissociation)
                self._write_dict_block(f, 'Reactions',             idata.reactions)
                self._write_dict_block(f, 'Derived (computed)',    derived)
 
                f.write("## (2) User Selections\n\n")
                self._write_dict_block(f, 'Run configuration', user_selections)
 
                f.write("## (3) Solver Configuration\n\n")
                self._write_dict_block(f, 'Solver settings', sc_flat)
 
                f.write("## (4) Run Statistics\n\n")
                self._write_dict_block(f, 'Runtime and machine', run_stats)
            print(f"  ✓ provenance.md")
        except Exception as exc:
            print(f"  ✗ provenance.md: {type(exc).__name__}: {exc}")
 
        # ── 2. Binary results (.npy) — write FIRST so the raw data is safe ───
        # even if everything else fails.  These are always present in `results`
        # for any non-empty run.
        try:
            t_arr = results.get('t')
            y_arr = results.get('y')
            if t_arr is not None and len(t_arr) > 0:
                np.save(str(out_dir / 'results_t.npy'),  t_arr)
                if y_arr is not None:
                    np.save(str(out_dir / 'results_y.npy'),  y_arr)
                print(f"  ✓ results_t.npy / results_y.npy")
            else:
                print(f"  ✗ results_t.npy: no time points in results")
        except Exception as exc:
            print(f"  ✗ results_*.npy: {type(exc).__name__}: {exc}")
 
        # ── 3. Summary CSV ───────────────────────────────────────────────────
        try:
            import csv
            row = post_process.summary_csv_row(results, self.input_data,
                                               solver_label=f"{sm}/{po}")
            with open(out_dir / 'summary.csv', 'w', newline='', encoding='utf-8') as f:
                writer = csv.DictWriter(f, fieldnames=list(row.keys()))
                writer.writeheader()
                writer.writerow(row)
            print(f"  ✓ summary.csv")
        except Exception as exc:
            print(f"  ✗ summary.csv: {type(exc).__name__}: {exc}")
 
        # ── 4. Diagnostics — guard each field access individually ────────────
        def _last(key, fmt='{:.6e}', mul=1.0):
            arr = results.get(key)
            if arr is None or len(arr) == 0:
                return '<unavailable>'
            try:
                return fmt.format(float(arr[-1]) * mul)
            except Exception:
                return '<unavailable>'
 
        try:
            diag_path = out_dir / 'diagnostics.txt'
            with open(diag_path, 'w', encoding='utf-8') as f:
                f.write(f"# Diagnostics for {label}\n\n")
                f.write("## Rate constants\n")
                f.write(getattr(self, '_diag_text', '') + '\n\n')
                f.write("## Conservation diagnostics (final time step)\n")
                f.write(f"delta_FP = {_last('delta_FP')}  (Frenkel pair)\n")
                f.write(f"delta_He = {_last('delta_He')}  (He balance)\n\n")
                f.write("## Key results\n")
                f.write(f"Final dose:       {_last('dose', '{:.4e}')} dpa\n")
                f.write(f"Swelling (final): {_last('swelling', '{:.6f}', 100.0)} %\n")
                f.write(f"C_He_tot (final): {_last('C_He_tot', '{:.3e}')} m^-3\n")
                f.write(f"mean_n_i (final): {_last('mean_n_i', '{:.2f}')}\n")
                f.write(f"mean_n_v (final): {_last('mean_n_v', '{:.2f}')}\n\n")
                prog = getattr(self, '_progress_rows', None)
                if prog:
                    f.write(f"## Time-step diagnostics ({len(prog)} rows)\n")
                    keys = list(prog[0].keys())
                    f.write('\t'.join(keys) + '\n')
                    for row in prog:
                        f.write('\t'.join(f"{row.get(k, 0.0):.6e}" for k in keys) + '\n')
            print(f"  ✓ diagnostics.txt")
        except Exception as exc:
            print(f"  ✗ diagnostics.txt: {type(exc).__name__}: {exc}")
 
        # ── 5. Plots — most likely to fail on degenerate partial data ────────
        try:
            from . import visualization
            visualization.save_all_plots(results, self.input_data,
                                         str(plot_dir), label=f"{sm}/{po}",
                                         rate_eq_obj=self.rate_equations)
            print(f"  ✓ plots/")
        except Exception as exc:
            print(f"  ✗ plots/: {type(exc).__name__}: {exc}  "
                  f"(numerical artifacts above are still safe)")
 
        print(f"Output written to: {out_dir}")
        # Remembered so a caller that did not choose the name can still find
        # the directory (the verification runner attaches provenance.md and
        # summary.csv from here to the record it pushes to the campaign board).
        self._last_output_dir = out_dir
        return out_dir