gammalooprs/settings/

runtime.rs

use std::{
    collections::{BTreeMap, BTreeSet},
    fmt::Display,
};

use bincode_trait_derive::{Decode, Encode};
use eyre::Result;
use linnet::half_edge::involution::EdgeVec;
use schemars::JsonSchema;
use serde::{Deserialize, Deserializer, Serialize, Serializer};
use spenso::algebra::{algebraic_traits::IsZero, complex::Complex};
use tracing::warn;
use typed_index_collections::TiVec;

use crate::{
    DependentMomentaConstructor, GammaLoopContext,
    cff::esurface::Esurface,
    graph::LoopMomentumBasis,
    integrands::process::SamplingMapDefinition,
    integrands::process::evaluators::EvaluatorMethod,
    momentum::{Helicity, RotationMethod, sample::ExternalIndex, signature::SignatureLike},
    observables::ObservableFileFormat,
    settings::runtime::kinematic::{Externals, improvement::generate_default_momenta},
    utils::{
        ApproxEq, DEFAULT_ESURFACE_EXISTENCE_THRESHOLD, F, FloatLike, format_uncertainty,
        serde_utils::{
            _default_rotation_axis, _default_stability_levels, IsDefault,
            deserialize_nonnegative_finite_f64, is_default_esurface_existence_threshold,
            is_default_rotation_axis, is_default_stability_levels, is_false, is_float, is_true,
            is_u64, is_usize, show_defaults_helper,
        },
    },
};
use symbolica::domains::{float::Real, integer::Integer, rational::Rational};

use super::{RuntimeSettings, global::OrientationPattern};

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Default, Encode, Decode, JsonSchema)]
#[serde(default)]
pub struct RuntimeModelSettings {
    #[serde(flatten, skip_serializing_if = "IsDefault::is_default")]
    pub external_parameters: BTreeMap<String, (F<f64>, F<f64>)>,
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[serde(default, deny_unknown_fields)]
pub struct SubtractionSettings {
    /// Configuration for local ultraviolet counterterms in the numerical integrand.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub local_ct_settings: LocalCounterTermSettings,
    /// Configuration for analytically integrated counterterm contributions.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub integrated_ct_settings: IntegratedCounterTermSettings,
    /// Strategy used to choose centers shared by overlapping threshold surfaces.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub overlap_settings: OverlapSettings,
    /// Dimensionless tolerance used to compare the energy-squared E-surface invariant margin
    /// against `esurface_existence_threshold * E_cm^2` at runtime.
    #[serde(
        deserialize_with = "deserialize_nonnegative_finite_f64",
        skip_serializing_if = "is_default_esurface_existence_threshold"
    )]
    #[schemars(range(min = 0.0))]
    pub esurface_existence_threshold: f64,
    /// Multiplier of `epsilon(T) * E_cm` used as the direct residual limit for
    /// threshold-counterterm radial roots at numerical precision `T`. A residual-limited
    /// higher-precision root may instead be accepted against the corresponding lower-precision
    /// limit after satisfying the cross-precision consistency checks.
    #[serde(
        deserialize_with = "deserialize_nonnegative_finite_f64",
        skip_serializing_if = "is_float::<64>"
    )]
    #[schemars(range(min = 0.0))]
    pub radial_root_residual_tolerance: f64,
    /// Disable local threshold subtraction while leaving the original integrand unchanged.
    #[serde(skip_serializing_if = "is_false")]
    pub disable_threshold_subtraction: bool,
}

impl Default for SubtractionSettings {
    fn default() -> Self {
        Self {
            local_ct_settings: LocalCounterTermSettings::default(),
            integrated_ct_settings: IntegratedCounterTermSettings::default(),
            overlap_settings: OverlapSettings::default(),
            esurface_existence_threshold: DEFAULT_ESURFACE_EXISTENCE_THRESHOLD,
            radial_root_residual_tolerance: 64.0,
            disable_threshold_subtraction: false,
        }
    }
}

#[derive(Copy, Clone)]
pub struct LockedRuntimeSettings<'a>(&'a RuntimeSettings);
impl<'a> From<&'a RuntimeSettings> for LockedRuntimeSettings<'a> {
    fn from(value: &'a RuntimeSettings) -> Self {
        LockedRuntimeSettings(value)
    }
}

impl<'a> From<LockedRuntimeSettings<'a>> for RuntimeSettings {
    fn from(value: LockedRuntimeSettings) -> Self {
        value.0.clone()
    }
}

impl<'a> LockedRuntimeSettings<'a> {
    pub(crate) fn helicities(&self) -> &[Helicity] {
        self.0.kinematics.externals.get_helicities()
    }

    /// If no external momenta are set in the settings, generate a default kinematic point that satisfies the phase-space constraints
    pub(crate) fn into_with_modified_kinematics(
        self,
        external_signature: &SignatureLike<ExternalIndex>,
        external_masses: &TiVec<ExternalIndex, F<f64>>,
    ) -> Result<RuntimeSettings> {
        if external_signature.is_empty() {
            Ok(self.into())
        } else {
            match &self.0.kinematics.externals {
                Externals::Constant {
                    momenta,
                    helicities,
                    ..
                } => {
                    if momenta.is_empty() && helicities.is_empty() {
                        warn!(
                            "No external kinematics were provided; using default-generated external momenta and helicities."
                        );
                        let new_externals = generate_default_momenta(
                            external_masses,
                            external_signature,
                            &F(self.0.kinematics.e_cm),
                        )?;

                        let mut new_settigns = self.0.clone();
                        new_settigns.kinematics.externals = new_externals;
                        Ok(new_settigns)
                    } else {
                        Ok(self.into())
                    }
                }
            }
        }
    }

    pub(crate) fn existence_check(
        &self,
        esurface: &Esurface,
        masses: &EdgeVec<F<f64>>,
        external_signature: &SignatureLike<ExternalIndex>,
        lmb: &LoopMomentumBasis,
        esurface_existence_threshold: f64,
    ) -> bool {
        let dependent_momenta_constructor =
            DependentMomentaConstructor::Amplitude(external_signature);

        esurface
            .classify_existence(
                &self
                    .0
                    .kinematics
                    .externals
                    .get_dependent_externals(dependent_momenta_constructor)
                    .unwrap(),
                lmb,
                masses,
                &F(self.0.kinematics.e_cm),
                &F(esurface_existence_threshold),
            )
            .is_existing()
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object))]
#[derive(
    Debug, Copy, Clone, Deserialize, Serialize, Encode, Decode, PartialEq, Eq, JsonSchema, Default,
)]
#[serde(rename_all = "snake_case")]
pub enum IntegralUnit {
    #[default]
    Auto,
    Picobarn,
    Femtobarn,
    Attobarn,
    Millibarn,
    None,
}

impl IntegralUnit {
    pub(crate) fn resolve_for_cross_section(self, n_initial_state_particles: usize) -> Self {
        match self {
            Self::Auto if n_initial_state_particles > 1 => Self::Picobarn,
            Self::Auto => Self::None,
            explicit => explicit,
        }
    }

    pub(crate) fn relative_to_picobarn_factor<T: FloatLike>(self, one: &F<T>) -> Option<F<T>> {
        match self {
            Self::Picobarn => Some(one.one()),
            Self::Femtobarn => Some(one.from_i64(1_000)),
            Self::Attobarn => Some(one.from_i64(1_000_000)),
            Self::Millibarn => Some(one.one() / one.from_i64(1_000_000_000)),
            Self::None => None,
            Self::Auto => unreachable!("integral unit must be resolved before conversion"),
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Deserialize, Serialize, Encode, Decode, PartialEq, JsonSchema)]
#[trait_decode(trait= GammaLoopContext)]
#[serde(default, deny_unknown_fields)]
pub struct GeneralSettings {
    /// Evaluate through the loop-tree-duality representation where the selected integrand supports it.
    #[serde(skip_serializing_if = "is_false")]
    pub use_ltd: bool,
    /// Evaluator backend used for generated integrands.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub evaluator_method: EvaluatorMethod,
    /// Optional orientation pattern restricting which causal-flow orientations are evaluated.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub orientation_pat: OrientationPattern,
    /// Load a previously compiled causal-flow function instead of rebuilding it for this run.
    #[serde(skip_serializing_if = "is_false")]
    pub load_compiled_cff: bool,
    /// Cache reusable evaluation data between phase-space samples.
    #[serde(skip_serializing_if = "is_false")]
    pub enable_cache: bool,
    /// Emit additional cache consistency diagnostics while evaluating samples.
    #[serde(skip_serializing_if = "is_false")]
    pub debug_cache: bool,
    /// Ultraviolet reference mass used by local counterterms, in the model's energy units.
    #[serde(skip_serializing_if = "is_float::<1000>")]
    pub m_uv: f64,
    /// Localization scale applied to ultraviolet counterterms, in the model's energy units.
    #[serde(skip_serializing_if = "is_float::<1000>")]
    pub renormalization_localization_scale: f64,
    /// Renormalization scale used by running parameters and scale-dependent terms.
    #[serde(skip_serializing_if = "is_float::<1000>")]
    pub mu_r: f64,
    /// Shared energy scale for numerator reconstruction when enabled during generation.
    #[serde(skip_serializing_if = "is_float::<1>")]
    pub numerator_sampling_scale: f64,
    /// Inverse temperature used to evaluate thermal distribution functions.
    #[serde(skip_serializing_if = "is_float::<1>")]
    pub inverse_temperature: f64,
    /// Values assigned, in order, to extra symbolic parameters expected by the evaluator.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub additional_param_values: Vec<f64>,
    /// Unit used to report integrated cross-sections; `auto` selects picobarns for scattering.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub integral_unit: IntegralUnit,
    /// Omit the incoming-state flux factor from the returned integrand value.
    #[serde(skip_serializing_if = "is_false")]
    pub disable_flux_factor: bool,
    /// Return generated events in addition to accumulating integration statistics.
    #[serde(skip_serializing_if = "is_false")]
    pub generate_events: bool,
    /// Persist every configured additional weight on each generated event.
    #[serde(skip_serializing_if = "is_false")]
    pub store_additional_weights_in_event: bool,
}

impl Default for GeneralSettings {
    fn default() -> Self {
        Self {
            evaluator_method: EvaluatorMethod::default(),
            use_ltd: false,
            load_compiled_cff: false,
            enable_cache: false,
            debug_cache: false,
            orientation_pat: OrientationPattern::default(),
            m_uv: 1000.0,
            renormalization_localization_scale: 1000.0,
            mu_r: 1000.0,
            numerator_sampling_scale: 1.0,
            inverse_temperature: 1.0,
            additional_param_values: vec![],
            integral_unit: IntegralUnit::Auto,
            disable_flux_factor: false,
            generate_events: false,
            store_additional_weights_in_event: false,
        }
    }
}

impl GeneralSettings {
    pub(crate) fn mu_r_sq(&self) -> f64 {
        self.mu_r * self.mu_r
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object))]
#[derive(
    Debug, Copy, Clone, PartialEq, Deserialize, Default, Serialize, Encode, Decode, JsonSchema,
)]
// #[trait_decode(trait= GammaLoopContext)]
#[serde(deny_unknown_fields)]
pub enum IntegratedPhase {
    #[serde(rename = "real")]
    #[default]
    Real,
    #[serde(rename = "imag")]
    Imag,
    #[serde(rename = "both")]
    Both,
}

pub mod kinematic;

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Deserialize, Serialize, Encode, Decode, PartialEq, JsonSchema)]
#[serde(default, deny_unknown_fields)]
pub struct ObservablesOutputSettings {
    /// Output formats written for configured observables; duplicates and `none` are ignored.
    #[serde(
        default = "default_observable_output_formats",
        skip_serializing_if = "is_default_observable_output_formats"
    )]
    pub format: Vec<ObservableFileFormat>,
}

fn default_observable_output_formats() -> Vec<ObservableFileFormat> {
    vec![ObservableFileFormat::Json]
}

fn is_default_observable_output_formats(formats: &Vec<ObservableFileFormat>) -> bool {
    show_defaults_helper(formats.as_slice() == default_observable_output_formats().as_slice())
}

impl Default for ObservablesOutputSettings {
    fn default() -> Self {
        Self {
            format: default_observable_output_formats(),
        }
    }
}

impl ObservablesOutputSettings {
    pub(crate) fn resolved_formats(&self) -> Vec<ObservableFileFormat> {
        let mut resolved_formats = Vec::new();
        for format in &self.format {
            if *format == ObservableFileFormat::None || resolved_formats.contains(format) {
                continue;
            }
            resolved_formats.push(*format);
        }
        resolved_formats
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Deserialize, Serialize, Encode, Decode, PartialEq, JsonSchema)]
#[serde(default, deny_unknown_fields)]
pub struct IntegratorSettings {
    /// Number of bins maintained along each continuous integration dimension.
    #[serde(skip_serializing_if = "is_usize::<64>")]
    pub n_bins: usize,
    /// Optional per-iteration bin counts, overriding the fixed `n_bins` value as training proceeds.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub bin_number_evolution: Option<Vec<usize>>,
    /// Minimum accumulated samples before an adaptive grid or discrete probability is updated.
    #[serde(skip_serializing_if = "is_usize::<1000>")]
    pub min_samples_for_update: usize,
    /// Number of samples in the first Monte Carlo iteration.
    #[serde(skip_serializing_if = "is_usize::<100000>")]
    pub n_start: usize,
    /// Samples added to each successive Monte Carlo iteration.
    #[serde(skip_serializing_if = "is_usize::<10000>")]
    pub n_increase: usize,
    /// Hard upper bound on the total number of integrand evaluations.
    #[serde(skip_serializing_if = "is_usize::<10000000000>")]
    pub n_max: usize,
    /// Stop once the estimated relative uncertainty reaches this value, if configured.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub target_relative_accuracy: Option<f64>,
    /// Stop once the estimated absolute uncertainty reaches this value, if configured.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub target_absolute_accuracy: Option<f64>,
    /// Complex component to integrate: real, imaginary, or both.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub integrated_phase: IntegratedPhase,
    /// Update rate for learned probabilities over discrete sampling choices.
    #[serde(skip_serializing_if = "is_float::<1>")]
    pub discrete_dim_learning_rate: f64,
    /// Update rate for adaptive grids over continuous integration dimensions.
    #[serde(skip_serializing_if = "is_float::<1>")]
    pub continuous_dim_learning_rate: f64,
    /// Train adaptive distributions on iteration averages instead of individual sample weights.
    #[serde(skip_serializing_if = "is_false")]
    pub train_on_avg: bool,
    /// Include maximum-weight diagnostics in integration progress and result data.
    #[serde(skip_serializing_if = "is_true")]
    pub show_max_wgt_info: bool,
    /// Cap the ratio used when turning sampled weights into adaptive probabilities.
    #[serde(skip_serializing_if = "is_float::<30>")]
    pub max_prob_ratio: f64,
    /// Seed for deterministic initialization of Monte Carlo random-number streams.
    #[serde(skip_serializing_if = "is_u64::<69>")]
    pub seed: u64,
    /// File formats used when writing observable histograms and distributions.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub observables_output: ObservablesOutputSettings,
}

impl Default for IntegratorSettings {
    fn default() -> Self {
        Self {
            n_bins: 64,
            bin_number_evolution: None,
            min_samples_for_update: 1000,
            n_start: 100000,
            n_increase: 10000,
            n_max: 10000000000,
            target_relative_accuracy: None,
            target_absolute_accuracy: None,
            integrated_phase: IntegratedPhase::default(),
            discrete_dim_learning_rate: 1.0,
            continuous_dim_learning_rate: 1.0,
            train_on_avg: false,
            show_max_wgt_info: true,
            max_prob_ratio: 30.0,
            seed: 69,
            observables_output: ObservablesOutputSettings::default(),
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Deserialize, Serialize, PartialEq, Encode, Decode, JsonSchema)]
#[serde(default, deny_unknown_fields)]
pub struct ParameterizationSettings {
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub mode: ParameterizationMode,
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub mapping: ParameterizationMapping,
    #[serde(skip_serializing_if = "is_float::<1>")]
    pub b: f64,
    #[serde(skip_serializing_if = "is_float::<1>")]
    pub power: f64,
    #[serde(skip_serializing_if = "BTreeMap::is_empty")]
    pub lmb_basis_ids: BTreeMap<String, Vec<usize>>,
    /// Graph-aware channels consumed by the production sampler. Keeping their
    /// selection with the parameterization preserves the canonical catalogue
    /// through the `SamplingSettings` enum round trip.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub sampling_channels: SamplingChannelSelection,
}

impl Default for ParameterizationSettings {
    fn default() -> Self {
        Self {
            b: 1.0,
            power: 1.0,
            mode: ParameterizationMode::default(),
            mapping: ParameterizationMapping::default(),
            lmb_basis_ids: BTreeMap::new(),
            sampling_channels: SamplingChannelSelection::default(),
        }
    }
}

#[derive(Serialize, Deserialize, Default, Debug, Clone)]
#[serde(default, deny_unknown_fields)]
pub struct IntegralEstimate {
    pub neval: usize,
    pub real_zero: usize,
    pub im_zero: usize,
    pub result: Complex<F<f64>>,
    pub error: Complex<F<f64>>,
    pub real_chisq: F<f64>,
    pub im_chisq: F<f64>,
}

impl Display for IntegralEstimate {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        // Format real and imaginary parts with proper uncertainty formatting
        let real_formatted = if self.error.re.0.abs() > 0.0 && f.alternate() {
            format_uncertainty(self.result.re, self.error.re)
        } else {
            // If no error, use formatter precision or default to scientific notation
            match f.precision() {
                Some(p) => format!("{:.prec$e}", self.result.re.0, prec = p),
                None => format!("{:e}", self.result.re.0),
            }
        };

        let im_formatted = if self.error.im.0.abs() > 0.0 && f.alternate() {
            format_uncertainty(self.result.im, self.error.im)
        } else {
            // If no error, use formatter precision or default to scientific notation
            match f.precision() {
                Some(p) => format!("{:.prec$e}", self.result.im.0, prec = p),
                None => format!("{:e}", self.result.im.0),
            }
        };

        // Display the result in a readable format
        if self.result.im.0.abs() > 0.0 {
            if f.alternate() {
                write!(
                    f,
                    "{} + {}i (neval: {}, real_chisq: {:.3}, im_chisq: {:.3})",
                    real_formatted, im_formatted, self.neval, self.real_chisq.0, self.im_chisq.0
                )
            } else {
                write!(f, "{} + {}i", real_formatted, im_formatted,)
            }
        } else if f.alternate() {
            write!(
                f,
                "{} (neval: {}, real_chisq: {:.3})",
                real_formatted, self.neval, self.real_chisq.0
            )
        } else {
            write!(f, "{}", real_formatted)
        }
    }
}

impl IntegralEstimate {
    pub fn is_compatible_with_target(&self, target: Complex<F<f64>>, sigma: i8) -> bool {
        let res_norm = self.result.norm().re;
        let tolerance = if res_norm.is_zero() {
            self.error.norm().re
        } else {
            (self.error).norm().re * F(sigma as f64) / res_norm
        };

        if tolerance.0 > 0.1 {
            warn!("Tolerance larger than 10%! {}", tolerance)
        }
        // println!("Tol: {}", tolerance.0);
        self.result.approx_eq(&target, &tolerance)
    }

    pub fn is_compatible_with_result(&self, other: &Self, sigma: i8) -> bool {
        let combined_error = Complex::new(
            (self.error.re * self.error.re + other.error.re * other.error.re).sqrt(),
            (self.error.im * self.error.im + other.error.im * other.error.im).sqrt(),
        );
        let delta = self.result - other.result;
        let delta_norm = delta.norm().re;
        let tolerance = if delta_norm.is_zero() {
            combined_error.norm().re
        } else {
            combined_error.norm().re * F(sigma as f64) / delta_norm
        };

        if tolerance.0 > 0.1 {
            warn!("Tolerance larger than 10%! {}", tolerance)
        }

        self.result.approx_eq(&other.result, &tolerance)
    }
}

#[derive(Serialize, Deserialize, Default, Debug, Clone)]
#[serde(default, deny_unknown_fields)]
pub struct IntegrationTableComponentResult {
    pub component: String,
    pub value: F<f64>,
    pub error: F<f64>,
    pub relative_error_percent: Option<f64>,
    pub chi_sq_per_dof: f64,
    pub target_delta_sigma: Option<f64>,
    pub target_delta_percent: Option<f64>,
    pub max_weight_impact: f64,
}

#[derive(Serialize, Deserialize, Default, Debug, Clone)]
#[serde(default, deny_unknown_fields)]
pub struct IntegrationStatisticsSnapshot {
    pub num_evals: usize,
    pub average_total_time_seconds: f64,
    pub average_parameterization_time_seconds: f64,
    pub average_integrand_time_seconds: f64,
    pub average_evaluator_time_seconds: f64,
    pub average_observable_time_seconds: f64,
    pub average_integrator_time_seconds: f64,
    pub f64_percentage: f64,
    pub f128_percentage: f64,
    pub arb_percentage: f64,
    pub nan_percentage: f64,
    pub nan_or_unstable_percentage: f64,
    pub generated_event_count: usize,
    pub accepted_event_count: usize,
    pub selection_efficiency_percentage: Option<f64>,
}

#[derive(Serialize, Deserialize, Default, Debug, Clone)]
#[serde(default, deny_unknown_fields)]
pub struct MaxWeightInfoEntry {
    pub component: String,
    pub sign: String,
    pub max_eval: F<f64>,
    pub coordinates: Option<String>,
}

#[derive(Serialize, Deserialize, Default, Debug, Clone)]
#[serde(default, deny_unknown_fields)]
pub struct DiscreteCoordinate {
    pub axis_label: String,
    pub bin_index: usize,
    pub bin_label: Option<String>,
}

#[derive(Serialize, Deserialize, Default, Debug, Clone)]
#[serde(default, deny_unknown_fields)]
pub struct DiscreteBreakdownEntry {
    pub bin_index: usize,
    pub bin_label: Option<String>,
    pub pdf: F<f64>,
    pub value: F<f64>,
    pub error: F<f64>,
    pub chi_sq: F<f64>,
    pub processed_samples: usize,
}

#[derive(Serialize, Deserialize, Default, Debug, Clone)]
#[serde(default, deny_unknown_fields)]
pub struct DiscreteBreakdown {
    pub axis_label: String,
    pub fixed_coordinates: Vec<DiscreteCoordinate>,
    pub entries: Vec<DiscreteBreakdownEntry>,
}

#[derive(Serialize, Deserialize, Default, Debug, Clone)]
#[serde(default, deny_unknown_fields)]
pub struct ComponentDiscreteBreakdown {
    pub re: Option<DiscreteBreakdown>,
    pub im: Option<DiscreteBreakdown>,
}

#[derive(Serialize, Deserialize, Default, Debug, Clone)]
#[serde(default, deny_unknown_fields)]
pub struct AbsoluteIntegrationResult {
    pub integral: IntegralEstimate,
    pub table_results: Vec<IntegrationTableComponentResult>,
    pub max_weight_info: Vec<MaxWeightInfoEntry>,
    pub grid_breakdown: ComponentDiscreteBreakdown,
}

#[derive(Serialize, Deserialize, Default, Debug, Clone)]
#[serde(default, deny_unknown_fields)]
pub struct SlotIntegrationResult {
    pub key: String,
    pub process: String,
    pub integrand: String,
    pub target: Option<Complex<F<f64>>>,
    pub integral: IntegralEstimate,
    pub table_results: Vec<IntegrationTableComponentResult>,
    pub integration_statistics: IntegrationStatisticsSnapshot,
    pub max_weight_info: Vec<MaxWeightInfoEntry>,
    pub grid_breakdown: ComponentDiscreteBreakdown,
    pub absolute: AbsoluteIntegrationResult,
}

#[derive(Serialize, Deserialize, Default, Debug, Clone)]
#[serde(default, deny_unknown_fields)]
pub struct IntegrationResult {
    pub slots: Vec<SlotIntegrationResult>,
}

impl IntegrationResult {
    pub fn slot(&self, key: &str) -> Option<&SlotIntegrationResult> {
        self.slots.iter().find(|slot| slot.key == key)
    }

    pub fn single_slot(&self) -> Option<&SlotIntegrationResult> {
        (self.slots.len() == 1).then(|| &self.slots[0])
    }

    pub fn single_slot_integral(&self) -> Option<&IntegralEstimate> {
        self.single_slot().map(|slot| &slot.integral)
    }

    pub fn is_compatible_with_target(&self, target: Complex<F<f64>>, n_sigma: i8) -> bool {
        self.single_slot_integral()
            .is_some_and(|integral| integral.is_compatible_with_target(target, n_sigma))
    }
}

impl Display for IntegrationResult {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        if let Some(slot) = self.single_slot() {
            return slot.integral.fmt(f);
        }

        for (index, slot) in self.slots.iter().enumerate() {
            if index > 0 {
                writeln!(f)?;
            }
            write!(f, "{}: {}", slot.key, slot.integral)?;
        }
        Ok(())
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::utils::serde_utils::ShowDefaultsGuard;
    use spenso::algebra::complex::Complex;

    #[test]
    fn sampling_lu_h_profile_shorthand_and_table_share_one_definition() {
        let shorthand: SamplingChannelDefinition =
            toml::from_str("around='phase_space(cut(2,3))'\nradial_profile='lu_h'").unwrap();
        let table: SamplingChannelDefinition = toml::from_str("around='phase_space(cut(2,3))'\nradial_profile={kind='lu_h',approximation='log_logistic',broad_fraction=0.02}").unwrap();
        assert_eq!(shorthand, table);
        assert_eq!(
            shorthand.radial_profile,
            Some(SamplingRadialProfile::default())
        );
        let exported = toml::to_string(&shorthand).unwrap();
        assert_eq!(
            toml::from_str::<SamplingChannelDefinition>(&exported).unwrap(),
            shorthand
        );
        for invalid in [
            "{kind='lu_h',broad_faction=0.1}",
            "'unknown'",
            "{kind='lu_h',approximation='exact'}",
        ] {
            assert!(
                toml::from_str::<SamplingChannelDefinition>(&format!("radial_profile={invalid}"))
                    .is_err()
            );
        }
        let _guard = ShowDefaultsGuard::new(true);
        let json = serde_json::to_value(SamplingRadialProfile::default()).unwrap();
        assert!(json.get("scale").is_some());
        assert!(json.get("shape").is_some());
    }

    #[test]
    fn h_function_normalizations_match_bessel_integrals_in_native_precision() {
        use symbolica::{
            atom::{Atom, AtomView},
            coefficient::CoefficientView,
            domains::float::{Float, RealLike},
            function,
            transcendental::bessel_k,
        };

        use crate::utils::{ArbPrec, SamplingFloat, f128};

        const ORACLE_PRECISION: u32 = 1200;
        // Symbolica's own arbitrary-precision evaluation of e^2 K_nu(2).
        let scaled_bessel_k = |twice_order: usize| {
            let value = function!(
                bessel_k(),
                Atom::num((twice_order as i64, 2)),
                Atom::num(Float::with_val(ORACLE_PRECISION, 2))
            );
            let AtomView::Num(number) = value.as_view() else {
                panic!("bessel_k did not evaluate: {value}");
            };
            let CoefficientView::Float(real, _) = number.get_coeff_view() else {
                panic!("bessel_k did not evaluate to a float: {value}");
            };
            real.to_float() * Float::with_val(ORACLE_PRECISION, 2).exp()
        };
        // The previously tabulated binary64 integrals anchor the Bessel orders.
        let sqrt_pi_half = std::f64::consts::PI.sqrt() / 2.0;
        let poly = [
            (None, sqrt_pi_half),
            (Some(1), 0.841_568_215_070_771_4),
            (Some(3), 1.033_476_847_068_688_6),
            (Some(4), 1.329_340_388_179_137),
            (Some(6), 2.880_237_507_721_463_7),
            (Some(7), 4.783_566_971_347_609),
            (Some(9), 16.225_745_976_182_285),
            (Some(10), 32.735_007_058_911_25),
            (Some(12), 155.837_465_922_583_42),
            (Some(13), 364.658_500_356_566_04),
            (Some(15), 2_257.637_553_015_473),
            (Some(16), 5_939.804_418_537_864),
        ];
        let left_right = [
            (Some(0), 2.066_953_694_137_377),
            (Some(1), 1.683_136_430_141_542_8),
            (Some(3), 3.750_090_124_278_92),
            (Some(4), 9.567_133_942_695_218),
            (Some(6), 139.373_101_752_153_5),
            (Some(7), 729.317_000_713_132_1),
            (Some(9), 32_336.242_742_929_753),
            (Some(10), 263_205.217_049_469),
            (Some(12), 2.427_503_717_893_097_5e7),
            (Some(13), 2.694_265_921_644_289e8),
            (Some(15), 4.261_555_045_314_143e10),
            (Some(16), 5.998_751_004_871_322e11),
        ];
        let tabulated = std::iter::once((HFunction::Exponential, None, sqrt_pi_half))
            .chain(poly.map(|(power, value)| (HFunction::PolyExponential, power, value)))
            .chain(
                left_right
                    .map(|(power, value)| (HFunction::PolyLeftRightExponential, power, value)),
            );
        fn check<T: FloatLike>(sigma: F<T>, settings: &HFunctionSettings, expected: &Float) {
            // Repeated calls must reuse only the unit integral, not a previous sigma.
            for multiplier in [1, 2, 3] {
                let scaled_sigma = &sigma * sigma.from_usize(multiplier);
                let expected =
                    expected.clone() * Float::with_val(ORACLE_PRECISION, multiplier as u32);
                let actual: Float = settings.normalization(&scaled_sigma).0.into();
                let relative = ((actual - &expected) / &expected).norm();
                // The K_0 series cancels about four bits before the upward recurrence.
                let epsilon: Float = sigma.epsilon().0.into();
                assert!(
                    relative <= epsilon * 64,
                    "{settings:?} in {}: relative error {relative}",
                    std::any::type_name::<T>()
                );
            }
        }
        for (function, power, binary64) in tabulated {
            let settings = HFunctionSettings {
                function: function.clone(),
                power,
                ..Default::default()
            };
            let p = power.unwrap_or(0);
            let expected = match function {
                HFunction::Exponential => scaled_bessel_k(1),
                HFunction::PolyExponential => scaled_bessel_k(p.abs_diff(1)),
                _ => scaled_bessel_k(2 * p.abs_diff(1)) * 2,
            };
            assert!(
                ((expected.to_f64() - binary64) / binary64).abs() < 4.0 * f64::EPSILON,
                "{settings:?}: {expected} != tabulated {binary64}"
            );
            // A sigma that is exact in binary keeps every precision comparable.
            let expected = expected * Float::with_val(ORACLE_PRECISION, 1.25);
            check(F(1.25_f64), &settings, &expected);
            check(F(f128::from_f64(1.25)), &settings, &expected);
            check(F(SamplingFloat::from_f64(1.25)), &settings, &expected);
            check(F(ArbPrec::from_f64(1.25)), &settings, &expected);
        }
        // The default profile keeps its exact sqrt(pi)*sigma/2 normalization.
        let sigma = F(ArbPrec::from_f64(1.25));
        assert_eq!(
            HFunctionSettings::default().normalization(&sigma),
            sigma.PI().sqrt() * &sigma / sigma.from_usize(2)
        );
    }

    #[test]
    fn test_integration_result_display() {
        let result = IntegralEstimate {
            neval: 1000000,
            real_zero: 0,
            im_zero: 0,
            result: Complex::new(F(1.23456789e-3), F(4.56789e-5)),
            error: Complex::new(F(1.2e-5), F(3.4e-7)),
            real_chisq: F(1.05),
            im_chisq: F(0.98),
        };

        let display_str = format!("{}", result);
        println!("Integration result display: {}", display_str);

        // Test with only real part
        let real_only_result = IntegralEstimate {
            neval: 500000,
            real_zero: 0,
            im_zero: 0,
            result: Complex::new(F(2.345e-2), F(0.0)),
            error: Complex::new(F(1.5e-4), F(0.0)),
            real_chisq: F(1.12),
            im_chisq: F(0.0),
        };

        let real_display_str = format!("{}", real_only_result);
        println!("Real-only result display: {}", real_display_str);

        // Test with custom precision
        let precision_str = format!("{:.2}", result);
        println!("With precision 2: {}", precision_str);
    }

    #[test]
    fn integrator_settings_target_accuracy_defaults_stay_optional() {
        let settings = IntegratorSettings::default();

        let hidden_defaults = toml::to_string(&settings).expect("serialize integrator defaults");
        assert!(!hidden_defaults.contains("target_relative_accuracy"));
        assert!(!hidden_defaults.contains("target_absolute_accuracy"));

        let _guard = ShowDefaultsGuard::new(true);
        let shown_defaults = toml::to_string(&settings).expect("serialize integrator defaults");
        assert!(shown_defaults.contains("n_bins"));
        assert!(!shown_defaults.contains("target_relative_accuracy"));
        assert!(!shown_defaults.contains("target_absolute_accuracy"));
    }

    #[test]
    fn integrator_settings_target_accuracy_serializes_when_configured() {
        let settings = IntegratorSettings {
            target_relative_accuracy: Some(0.05),
            target_absolute_accuracy: Some(1.0e-6),
            ..IntegratorSettings::default()
        };

        let serialized = toml::to_string(&settings).expect("serialize configured integrator");
        assert!(serialized.contains("target_relative_accuracy = 0.05"));
        assert!(serialized.contains("target_absolute_accuracy = "));

        let deserialized: IntegratorSettings =
            toml::from_str(&serialized).expect("deserialize configured integrator");
        assert_eq!(settings, deserialized);
    }

    #[test]
    fn stability_ecm_relative_tolerances_default_and_roundtrip() {
        let _guard = ShowDefaultsGuard::new(false);
        for settings in [
            StabilityLevelSetting::default_double(),
            StabilityLevelSetting::default_quad(),
            StabilityLevelSetting::default_arb(),
        ] {
            let serialized = toml::to_string(&settings).unwrap();
            assert!(!serialized.contains("ecm_relative_tolerance"));
            assert_eq!(
                toml::from_str::<StabilityLevelSetting>(&serialized).unwrap(),
                settings
            );
        }

        let settings = StabilityLevelSetting {
            ecm_relative_tolerance_for_re: 1e-100,
            ecm_relative_tolerance_for_im: 2e-90,
            ..StabilityLevelSetting::default_arb()
        };
        let serialized = toml::to_string(&settings).unwrap();
        assert_eq!(
            toml::from_str::<StabilityLevelSetting>(&serialized).unwrap(),
            settings
        );
        let schema = serde_json::to_value(schemars::schema_for!(StabilityLevelSetting)).unwrap();
        for field in [
            "ecm_relative_tolerance_for_re",
            "ecm_relative_tolerance_for_im",
        ] {
            assert_eq!(schema["properties"][field]["minimum"], 0.0);
        }

        let defaults: StabilitySettings = toml::from_str("").unwrap();
        assert_eq!(defaults.integrated_energy_dimension, None);
        assert!(
            !toml::to_string(&defaults)
                .unwrap()
                .contains("integrated_energy_dimension")
        );
        let declared = StabilitySettings {
            integrated_energy_dimension: Some(-6),
            levels: vec![settings],
            ..defaults
        };
        assert_eq!(
            toml::from_str::<StabilitySettings>(&toml::to_string(&declared).unwrap()).unwrap(),
            declared
        );
        drop(_guard);
        let _guard = ShowDefaultsGuard::new(true);
        for settings in [
            StabilityLevelSetting::default_double(),
            StabilityLevelSetting::default_quad(),
            StabilityLevelSetting::default_arb(),
        ] {
            let shown = serde_json::to_value(settings).unwrap();
            assert_eq!(shown["ecm_relative_tolerance_for_re"], 0.0);
            assert_eq!(shown["ecm_relative_tolerance_for_im"], 0.0);
        }
        assert_eq!(
            serde_json::to_value(&declared).unwrap()["integrated_energy_dimension"],
            -6
        );
    }

    #[test]
    fn stability_ecm_relative_tolerances_reject_invalid_values() {
        let _guard = ShowDefaultsGuard::new(false);
        let defaults = toml::to_string(&StabilityLevelSetting::default_arb()).unwrap();
        for field in [
            "ecm_relative_tolerance_for_re",
            "ecm_relative_tolerance_for_im",
        ] {
            for invalid in ["-1.0", "nan", "inf", "-inf"] {
                let text = format!("{defaults}\n{field} = {invalid}\n");
                assert!(
                    toml::from_str::<StabilityLevelSetting>(&text).is_err(),
                    "{field} = {invalid}"
                );
            }
        }
    }

    #[test]
    fn sampling_ose_and_named_weight_roundtrip() {
        let _guard = ShowDefaultsGuard::new(false);
        for (graphs, channels) in [
            ("summed", "summed"),
            ("monte_carlo", "summed"),
            ("monte_carlo", "monte_carlo"),
        ] {
            for (global, local) in [("ose", "map_density"), ("map_density", "ose")] {
                let text = format!(
                    r#"graphs = "{graphs}"
sampling_multichanneling = true
sampling_channels = "{channels}"
sampling_channel_weight = "{global}"
alpha = 2.5
default_channel_selection = ["auto:optimized_lmb"]
[channel_definitions.G.soft]
around = "lmb(1,2)"
parent_lmb = [1,2]
channel_weight = "{local}"
"#
                );
                let settings: SamplingSettings = toml::from_str(&text).unwrap();
                let selection = settings
                    .get_parameterization_settings()
                    .unwrap()
                    .sampling_channels;
                assert_eq!(selection.alpha, 2.5);
                assert_eq!(
                    selection.weight == SamplingChannelWeight::Ose,
                    global == "ose"
                );
                assert_eq!(
                    selection.channel_definitions["G"]["soft"].channel_weight
                        == Some(SamplingChannelWeight::Ose),
                    local == "ose"
                );
                let serialized = toml::to_string(&settings).unwrap();
                assert_eq!(
                    toml::from_str::<SamplingSettings>(&serialized).unwrap(),
                    settings
                );
            }
        }
        let settings: SamplingSettings =
            toml::from_str("sampling_multichanneling=true\nsampling_channel_weight='ose'").unwrap();
        let selection = settings
            .get_parameterization_settings()
            .unwrap()
            .sampling_channels;
        assert_eq!(selection.weight, SamplingChannelWeight::Ose);
        assert_eq!(selection.alpha, 3.0);
        assert_eq!(
            selection.default_channel_selection,
            vec!["auto:optimized_lmb"]
        );
        assert!(!toml::to_string(&settings).unwrap().contains("alpha"));
        drop(_guard);
        let _guard = ShowDefaultsGuard::new(true);
        assert_eq!(serde_json::to_value(settings).unwrap()["alpha"], 3.0);
    }

    #[test]
    fn sampling_ose_alpha_validation_and_zero_control() {
        for invalid in ["-1.0", "nan", "inf", "-inf"] {
            let text = format!("alpha={invalid}");
            assert!(toml::from_str::<SamplingSettings>(&text).is_err());
            assert!(toml::from_str::<SamplingChannelSelection>(&text).is_err());
        }
        for alpha in [-1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
            assert!(
                SamplingSettings::from_parser(SamplingSettingsParser {
                    alpha,
                    ..Default::default()
                })
                .is_err()
            );
        }
        let settings: SamplingSettings = toml::from_str(
            "sampling_multichanneling=true\nsampling_channel_weight='ose'\nalpha=0.0",
        )
        .unwrap();
        assert_eq!(
            settings
                .get_parameterization_settings()
                .unwrap()
                .sampling_channels
                .alpha,
            0.0
        );
        assert_eq!(
            toml::from_str::<SamplingSettings>(&toml::to_string(&settings).unwrap()).unwrap(),
            settings
        );
    }

    #[test]
    fn uv_localisation_smooth_sliver_defaults_and_roundtrips() {
        let defaults: UVLocalisationSettings = toml::from_str("").unwrap();
        assert!(!defaults.smooth_sliver);
        assert_eq!(defaults, UVLocalisationSettings::default());
        {
            let _guard = ShowDefaultsGuard::new(false);
            assert!(
                !toml::to_string(&defaults)
                    .unwrap()
                    .contains("smooth_sliver")
            );
            let settings = UVLocalisationSettings {
                smooth_sliver: true,
                ..defaults.clone()
            };
            let serialized = toml::to_string(&settings).unwrap();
            assert!(serialized.contains("smooth_sliver = true"));
            assert_eq!(
                toml::from_str::<UVLocalisationSettings>(&serialized).unwrap(),
                settings
            );
        }
        let _guard = ShowDefaultsGuard::new(true);
        assert!(
            toml::to_string(&defaults)
                .unwrap()
                .contains("smooth_sliver = false")
        );
    }

    #[test]
    fn integral_unit_auto_resolves_by_initial_state_count() {
        assert_eq!(
            IntegralUnit::Auto.resolve_for_cross_section(1),
            IntegralUnit::None
        );
        assert_eq!(
            IntegralUnit::Auto.resolve_for_cross_section(2),
            IntegralUnit::Picobarn
        );
    }

    #[test]
    fn integral_unit_relative_scalings_match_requested_barn_unit() {
        let one = F(1.0f64);
        assert_eq!(
            IntegralUnit::Picobarn.relative_to_picobarn_factor(&one),
            Some(F(1.0))
        );
        assert_eq!(
            IntegralUnit::Femtobarn.relative_to_picobarn_factor(&one),
            Some(F(1_000.0))
        );
        assert_eq!(
            IntegralUnit::Attobarn.relative_to_picobarn_factor(&one),
            Some(F(1_000_000.0))
        );
        assert_eq!(
            IntegralUnit::Millibarn.relative_to_picobarn_factor(&one),
            Some(F(1.0e-9))
        );
        assert_eq!(IntegralUnit::None.relative_to_picobarn_factor(&one), None);
    }

    #[test]
    fn general_settings_store_mu_r_and_compute_mu_r_sq_internally() {
        let settings = GeneralSettings {
            mu_r: 91.188,
            ..GeneralSettings::default()
        };

        assert_eq!(settings.mu_r, 91.188);
        assert!((settings.mu_r_sq() - 8315.251344).abs() < 1.0e-12);

        let serialized = toml::to_string(&settings).expect("serialize general settings");
        assert!(serialized.contains("mu_r = 91.188"));
        assert!(!serialized.contains("mu_r_sq"));
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Serialize, Deserialize, Debug, Clone, Encode, Decode, PartialEq, JsonSchema)]
#[serde(default, deny_unknown_fields)]
pub struct StabilitySettings {
    /// Energy dimension of the integrated physical quantity after flux
    /// normalization and before conversion to output units. Required when an
    /// E_cm-relative tolerance is enabled for a physical evaluation; reference
    /// functions carry their own known dimensions. Direct momentum evaluations
    /// subtract the dimension of their missing spatial integration measure.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub integrated_energy_dimension: Option<i32>,
    /// Rotations applied when comparing an unstable sample with an equivalent phase-space point.
    #[serde(skip_serializing_if = "is_default_rotation_axis")]
    pub rotation_axis: Vec<RotationSetting>,
    /// Ordered numerical-precision levels and tolerances used during stability escalation.
    #[serde(skip_serializing_if = "is_default_stability_levels")]
    pub levels: Vec<StabilityLevelSetting>,
    /// Compare rotated results after normalizing by the result magnitude.
    #[serde(skip_serializing_if = "is_true")]
    pub check_on_norm: bool,
    /// Escalate precision when a result is exactly zero at the current level.
    #[serde(skip_serializing_if = "is_false")]
    pub escalate_if_exact_zero: bool,
    /// Suppress non-final stability escalation for negligible weighted probes.
    /// When positive, a probe is allowed to remain at its current precision if
    /// its largest absolute weighted rotation result is below this fraction of
    /// an absolute-component integral estimate whose relative standard error
    /// is at most ten percent. Zero disables this waiver.
    #[serde(
        default,
        deserialize_with = "deserialize_nonnegative_finite_f64",
        skip_serializing_if = "is_float::<0>"
    )]
    #[schemars(range(min = 0.0))]
    pub min_abs_wgt_for_escalation: f64,
    /// Escalate when the loop-momentum norm exceeds this factor; a negative value disables it.
    #[serde(skip_serializing_if = "is_float::<-1>")]
    pub loop_momenta_norm_escalation_factor: f64,
    /// Optional diagnostics retained for rotated, escalated, and high-precision evaluations.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub recording: Option<StabilityRecordingSettings>,
}

impl Default for StabilitySettings {
    fn default() -> Self {
        Self {
            integrated_energy_dimension: None,
            rotation_axis: _default_rotation_axis(),
            levels: _default_stability_levels(),
            check_on_norm: true,
            escalate_if_exact_zero: false,
            min_abs_wgt_for_escalation: 0.0,
            loop_momenta_norm_escalation_factor: -1.0,
            recording: None,
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Serialize, Deserialize, Debug, Clone, Copy, Encode, Decode, PartialEq, JsonSchema)]
#[serde(default, deny_unknown_fields)]
pub struct StabilityRecordingSettings {
    /// Record the values obtained at rotated phase-space points.
    #[serde(skip_serializing_if = "is_false")]
    pub record_rotated_results: bool,
    /// Record the result produced at every attempted numerical-precision level.
    #[serde(skip_serializing_if = "is_false")]
    pub record_all_stability_levels: bool,
    /// Record evaluations escalated because of the loop-momentum norm threshold.
    #[serde(skip_serializing_if = "is_false")]
    pub record_loop_momenta_escalation: bool,
}

#[allow(clippy::derivable_impls)]
impl Default for StabilityRecordingSettings {
    fn default() -> Self {
        Self {
            record_rotated_results: false,
            record_all_stability_levels: false,
            record_loop_momenta_escalation: false,
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Serialize, Deserialize, Debug, Clone, Copy, Encode, Decode, PartialEq, JsonSchema)]
#[serde(deny_unknown_fields)]
pub struct StabilityLevelSetting {
    /// Numeric precision used for this escalation level.
    pub precision: Precision,
    /// Maximum accepted relative error in the real component before escalating.
    pub required_precision_for_re: f64,
    /// Maximum accepted relative error in the imaginary component before escalating.
    pub required_precision_for_im: f64,
    /// Dimensionless real-component probe disagreement relative to E_cm raised
    /// to the observable's declared energy dimension, with its output-unit
    /// conversion. Includes every sample factor and the outer grid weight;
    /// direct momentum densities account for the missing spatial measure.
    /// Zero disables this alternative to relative accuracy. This is not the
    /// integrator's target standard error.
    #[serde(
        default,
        deserialize_with = "deserialize_nonnegative_finite_f64",
        skip_serializing_if = "is_float::<0>"
    )]
    #[schemars(range(min = 0.0))]
    pub ecm_relative_tolerance_for_re: f64,
    /// Independent dimensionless imaginary-component tolerance, with the same
    /// E_cm-relative convention as `ecm_relative_tolerance_for_re`.
    #[serde(
        default,
        deserialize_with = "deserialize_nonnegative_finite_f64",
        skip_serializing_if = "is_float::<0>"
    )]
    #[schemars(range(min = 0.0))]
    pub ecm_relative_tolerance_for_im: f64,
    /// Ratio threshold for the large-weight escalation check; a non-positive value disables it.
    pub escalate_for_large_weight_threshold: f64,
}

impl StabilityLevelSetting {
    pub fn default_double() -> Self {
        Self {
            precision: Precision::Double,
            required_precision_for_re: 1e-5,
            required_precision_for_im: 1e-5,
            ecm_relative_tolerance_for_re: 0.0,
            ecm_relative_tolerance_for_im: 0.0,
            escalate_for_large_weight_threshold: 0.9,
        }
    }

    pub fn default_quad() -> Self {
        Self {
            precision: Precision::Quad,
            required_precision_for_re: 1e-5,
            required_precision_for_im: 1e-5,
            ecm_relative_tolerance_for_re: 0.0,
            ecm_relative_tolerance_for_im: 0.0,
            escalate_for_large_weight_threshold: -1.0,
        }
    }

    pub fn default_arb() -> Self {
        Self {
            precision: Precision::Arb,
            required_precision_for_re: 1e-5,
            required_precision_for_im: 1e-5,
            ecm_relative_tolerance_for_re: 0.0,
            ecm_relative_tolerance_for_im: 0.0,
            escalate_for_large_weight_threshold: -1.0,
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Serialize, Deserialize, Debug, Clone, Copy, PartialEq, Encode, Decode, JsonSchema)]
#[serde(tag = "type")]
#[serde(deny_unknown_fields)]
pub enum RotationSetting {
    #[serde(rename = "x")]
    Pi2X {},
    #[serde(rename = "y")]
    Pi2Y {},
    #[serde(rename = "z")]
    Pi2Z {},
    #[serde(rename = "none")]
    None {},
    #[serde(rename = "euler_angles")]
    EulerAngles {
        /// First Euler angle, in radians.
        alpha: f64,
        /// Second Euler angle, in radians.
        beta: f64,
        /// Third Euler angle, in radians.
        gamma: f64,
    },
}

impl Default for RotationSetting {
    fn default() -> Self {
        Self::Pi2Z {}
    }
}

impl RotationSetting {
    pub(crate) fn rotation_method(&self) -> RotationMethod {
        match self {
            Self::Pi2X {} => RotationMethod::Pi2X,
            Self::Pi2Y {} => RotationMethod::Pi2Y,
            Self::Pi2Z {} => RotationMethod::Pi2Z,
            Self::None {} => RotationMethod::Identity,
            Self::EulerAngles { alpha, beta, gamma } => {
                RotationMethod::EulerAngles(*alpha, *beta, *gamma)
            }
        }
    }

    pub(crate) fn _as_str(&self) -> String {
        match self {
            Self::Pi2X {} => "x".to_owned(),
            Self::Pi2Y {} => "y".to_owned(),
            Self::Pi2Z {} => "z".to_owned(),
            Self::None {} => "none".to_owned(),
            Self::EulerAngles { alpha, beta, gamma } => {
                format!("euler {} {} {}", alpha, beta, gamma)
            }
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object))]
#[derive(
    Serialize, Deserialize, Debug, Clone, PartialEq, Copy, Hash, Eq, Encode, Decode, JsonSchema,
)]
#[serde(deny_unknown_fields)]
pub enum Precision {
    Double,
    Quad,
    Arb,
}

impl Display for Precision {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Self::Double => write!(f, "f64"),
            Self::Quad => write!(f, "f128"),
            Self::Arb => write!(f, "arb"),
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Default, Encode, Decode, JsonSchema)]
#[serde(deny_unknown_fields)]
pub enum ParameterizationMode {
    #[serde(rename = "cartesian")]
    Cartesian,
    #[serde(rename = "spherical")]
    #[default]
    Spherical,
    #[serde(rename = "hyperspherical")]
    HyperSpherical,
    #[serde(rename = "hyperspherical_flat")]
    HyperSphericalFlat,
    #[serde(rename = "momentum_space")]
    MomentumSpace,
    #[serde(rename = "relative_spherical")]
    RelativeSpherical,
    #[serde(rename = "spherical_common_radial")]
    SphericalCommonRadial,
    #[serde(rename = "spherical_product_common_radial")]
    SphericalProductCommonRadial,
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Serialize, Deserialize, Default)]
#[serde(default, deny_unknown_fields)]
pub struct OutputMetadata {
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub model_name: String,
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub output_type: String,
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub contents: Vec<String>,
}

impl Display for ParameterizationMode {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            ParameterizationMode::Cartesian => write!(f, "cartesian"),
            ParameterizationMode::Spherical => write!(f, "spherical"),
            ParameterizationMode::RelativeSpherical => write!(f, "relative spherical"),
            ParameterizationMode::SphericalCommonRadial => {
                write!(f, "common-radial spherical")
            }
            ParameterizationMode::SphericalProductCommonRadial => {
                write!(f, "product/common-radial spherical")
            }
            ParameterizationMode::HyperSpherical => write!(f, "hyperspherical"),
            ParameterizationMode::HyperSphericalFlat => write!(f, "flat hyperspherical"),
            ParameterizationMode::MomentumSpace => write!(f, "momentum space"),
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object))]
#[derive(Debug, Clone, Deserialize, PartialEq, Default, Serialize, Encode, Decode, JsonSchema)]
#[serde(deny_unknown_fields)]
pub enum ParameterizationMapping {
    #[serde(rename = "log")]
    Log,
    #[serde(rename = "power")]
    Power,
    #[serde(rename = "linear")]
    #[default]
    Linear,
}

#[derive(Debug, Clone, PartialEq, Encode, Decode)]
#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
pub enum SamplingSettings {
    Default(ParameterizationSettings),
    MultiChanneling(MultiChannelingSettings),
    DiscreteGraphs(DiscreteGraphSamplingSettings),
}

#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[serde(default, deny_unknown_fields)]
pub struct SamplingSettingsParser {
    /// Whether graph contributions are summed deterministically or sampled as a discrete axis.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub graphs: SumMode,
    /// Optional graph-name allow-list for discrete graph sampling.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub graph_names: Vec<String>,
    /// Whether causal orientations are summed deterministically or sampled as a discrete axis.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub orientations: SumMode,
    /// Enable multichannel sampling over the resolved graph-aware channels.
    #[serde(skip_serializing_if = "is_false")]
    #[serde(rename = "sampling_multichanneling", alias = "lmb_multichanneling")]
    pub sampling_multichanneling: bool,
    /// Whether resolved sampling channels are summed or sampled discretely.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    #[serde(rename = "sampling_channels", alias = "lmb_channels")]
    pub sampling_channels: SumMode,
    /// Exponent of each selected loop edge's on-shell energy in an OSE channel score.
    #[serde(
        deserialize_with = "deserialize_nonnegative_finite_f64",
        skip_serializing_if = "is_float::<3>"
    )]
    #[schemars(range(min = 0.0))]
    pub alpha: f64,
    /// Rule used to construct the relative probability of each resolved sampling channel.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    #[serde(rename = "sampling_channel_weight", alias = "lmb_channel_weight")]
    pub sampling_channel_weight: SamplingChannelWeight,
    /// Coordinate system used to parameterize loop momenta.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub coordinate_system: CoordinateSystem,
    /// Mapping from the unit hypercube to unbounded radial coordinates.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub mapping: ParameterizationMapping,
    /// Scale parameter of the selected radial mapping.
    #[serde(skip_serializing_if = "is_float::<1>")]
    pub b: f64,
    /// Power parameter of the selected radial mapping.
    #[serde(skip_serializing_if = "is_float::<1>")]
    pub power: f64,
    /// Explicit loop-momentum-basis identifiers keyed by graph name.
    #[serde(skip_serializing_if = "BTreeMap::is_empty")]
    pub lmb_basis_ids: BTreeMap<String, Vec<usize>>,
    /// Default channel selectors for graphs without a graph-specific selection.
    #[serde(skip_serializing_if = "Vec::is_empty")]
    pub default_channel_selection: Vec<String>,
    /// Channel selectors keyed by graph name, overriding the default selection.
    #[serde(skip_serializing_if = "BTreeMap::is_empty")]
    pub channel_selection: BTreeMap<String, Vec<String>>,
    /// Named channel definitions keyed first by graph name, then by channel name.
    #[serde(skip_serializing_if = "BTreeMap::is_empty")]
    pub channel_definitions: BTreeMap<String, BTreeMap<String, SamplingChannelDefinition>>,
}

// Keep the old Python attribute spellings available while exposing the generic
// `sampling_*` names as the canonical fields.  Serde aliases above cover TOML
// and JSON cards; these accessors cover callers that construct parser settings
// directly through the Python bindings.
#[cfg(feature = "python")]
#[pyo3::pymethods]
impl SamplingSettingsParser {
    #[getter(lmb_multichanneling)]
    fn legacy_lmb_multichanneling(&self) -> bool {
        self.sampling_multichanneling
    }

    #[setter(lmb_multichanneling)]
    fn set_legacy_lmb_multichanneling(&mut self, value: bool) {
        self.sampling_multichanneling = value;
    }

    #[getter(lmb_channels)]
    fn legacy_lmb_channels(&self) -> SumMode {
        self.sampling_channels.clone()
    }

    #[setter(lmb_channels)]
    fn set_legacy_lmb_channels(&mut self, value: SumMode) {
        self.sampling_channels = value;
    }

    #[getter(lmb_channel_weight)]
    fn legacy_lmb_channel_weight(&self) -> SamplingChannelWeight {
        self.sampling_channel_weight
    }

    #[setter(lmb_channel_weight)]
    fn set_legacy_lmb_channel_weight(&mut self, value: SamplingChannelWeight) {
        self.sampling_channel_weight = value;
    }
}

impl Default for SamplingSettingsParser {
    fn default() -> Self {
        Self {
            graphs: SumMode::Summed,
            graph_names: Vec::new(),
            orientations: SumMode::Summed,
            sampling_multichanneling: false,
            sampling_channels: SumMode::Summed,
            alpha: 3.0,
            sampling_channel_weight: SamplingChannelWeight::default(),
            coordinate_system: CoordinateSystem::Spherical,
            mapping: ParameterizationMapping::Linear,
            b: 1.0,
            power: 1.0,
            lmb_basis_ids: BTreeMap::new(),
            default_channel_selection: Vec::new(),
            channel_selection: BTreeMap::new(),
            channel_definitions: BTreeMap::new(),
        }
    }
}

/// Weighting strategy for the resolved sampling-channel partition.
///
/// Exact inverse-map densities, on-shell-energy scores, and explicit positive
/// proxies share the same channel catalogue. `inverse_jacobian` selects the
/// same density score as `map_density`; `ose` is a distinct energy score.
#[derive(
    Debug, Default, Clone, Copy, Serialize, Deserialize, PartialEq, Eq, Encode, Decode, JsonSchema,
)]
#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
pub enum SamplingChannelWeight {
    #[serde(rename = "map_density")]
    #[default]
    MapDensity,
    #[serde(rename = "ose")]
    Ose,
    #[serde(rename = "singularity_proxy")]
    SingularityProxy,
    #[serde(rename = "inverse_jacobian")]
    InverseJacobian,
}

/// Radial proposal family for a physical-cut sampling channel.
#[derive(
    Debug, Default, Clone, Copy, Serialize, Deserialize, PartialEq, Eq, Encode, Decode, JsonSchema,
)]
#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[serde(rename_all = "snake_case")]
pub enum SamplingRadialProfileKind {
    #[default]
    LuH,
}

#[derive(
    Debug, Default, Clone, Copy, Serialize, Deserialize, PartialEq, Eq, Encode, Decode, JsonSchema,
)]
#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[serde(rename_all = "snake_case")]
pub enum SamplingRadialApproximation {
    #[default]
    LogLogistic,
}

/// A normalized proposal for the auxiliary LU scale. This changes sampling
/// only: the physical LU h-function and its raised-cut derivatives are intact.
#[derive(Debug, Clone, Serialize, PartialEq, Encode, Decode, JsonSchema)]
#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[schemars(with = "SamplingRadialProfileInput")]
pub struct SamplingRadialProfile {
    pub kind: SamplingRadialProfileKind,
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub approximation: SamplingRadialApproximation,
    pub broad_fraction: f64,
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub scale: Option<f64>,
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub shape: Option<f64>,
}

impl Default for SamplingRadialProfile {
    fn default() -> Self {
        Self {
            kind: SamplingRadialProfileKind::LuH,
            approximation: SamplingRadialApproximation::LogLogistic,
            broad_fraction: 0.02,
            scale: None,
            shape: None,
        }
    }
}

// Normalize the compact string and the explicit table to one runtime type.
// Export always retains the detailed form, including the proposal overrides.
#[derive(Deserialize, JsonSchema)]
#[serde(untagged, deny_unknown_fields)]
enum SamplingRadialProfileInput {
    Name(SamplingRadialProfileKind),
    Settings {
        kind: SamplingRadialProfileKind,
        #[serde(default)]
        approximation: SamplingRadialApproximation,
        #[serde(default = "SamplingRadialProfile::default_broad_fraction")]
        broad_fraction: f64,
        #[serde(default)]
        scale: Option<f64>,
        #[serde(default)]
        shape: Option<f64>,
    },
}

impl SamplingRadialProfile {
    fn default_broad_fraction() -> f64 {
        Self::default().broad_fraction
    }
}

impl<'de> Deserialize<'de> for SamplingRadialProfile {
    fn deserialize<D: serde::Deserializer<'de>>(deserializer: D) -> Result<Self, D::Error> {
        Ok(
            match SamplingRadialProfileInput::deserialize(deserializer)? {
                SamplingRadialProfileInput::Name(kind) => Self {
                    kind,
                    ..Self::default()
                },
                SamplingRadialProfileInput::Settings {
                    kind,
                    approximation,
                    broad_fraction,
                    scale,
                    shape,
                } => Self {
                    kind,
                    approximation,
                    broad_fraction,
                    scale,
                    shape,
                },
            },
        )
    }
}

/// A graph-scoped channel definition.  `parent_lmb` is deliberately explicit:
/// resolving a short expression against a graph is a later validation step,
/// while exported definitions must never lose the routing frame they use.
#[derive(Debug, Default, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[serde(default, deny_unknown_fields)]
pub struct SamplingChannelDefinition {
    pub around: String,
    /// Override the global channel score for this named channel. OSE scores
    /// require a complete ordinary LMB map; advanced maps keep their actual
    /// density unless an explicit compatible proxy is selected.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub channel_weight: Option<SamplingChannelWeight>,
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub radial_profile: Option<SamplingRadialProfile>,
    /// Ordered loop edges spanning the coordinates in which a surface or
    /// joint constraint is solved.  This is distinct from the edge set in
    /// `around`, which identifies the physical energy constraints.
    pub subspace_lmb: Vec<usize>,
    pub parent_lmb: Vec<usize>,
    pub on_cut: Vec<usize>,
    /// Optional positive singularity proxy used when the effective global or
    /// per-channel weight is `singularity_proxy`. The expression is
    /// parsed by Symbolica and evaluated in the complete master raw frame;
    /// its variables are `x0`, `x1`, ... in coordinate order.  A proxy is
    /// never inferred from the map, since doing so could hide an incomplete
    /// singularity model.
    #[serde(skip_serializing_if = "Option::is_none")]
    pub singularity_proxy: Option<String>,
}

/// Selection and user definitions for graph-aware sampling channels.
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[serde(default, deny_unknown_fields)]
pub struct SamplingChannelSelection {
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub weight: SamplingChannelWeight,
    /// Exponent of each selected LMB edge's on-shell energy in an OSE score.
    /// Zero gives a constant score; this is not a normalized proposal density.
    #[serde(
        deserialize_with = "deserialize_nonnegative_finite_f64",
        skip_serializing_if = "is_float::<3>"
    )]
    #[schemars(range(min = 0.0))]
    pub alpha: f64,
    pub default_channel_selection: Vec<String>,
    pub channel_selection: BTreeMap<String, Vec<String>>,
    pub channel_definitions: BTreeMap<String, BTreeMap<String, SamplingChannelDefinition>>,
}

impl Default for SamplingChannelSelection {
    fn default() -> Self {
        Self {
            weight: SamplingChannelWeight::default(),
            alpha: 3.0,
            default_channel_selection: Vec::new(),
            channel_selection: BTreeMap::new(),
            channel_definitions: BTreeMap::new(),
        }
    }
}

#[derive(Debug, Default, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
pub enum SumMode {
    #[serde(rename = "summed")]
    #[default]
    Summed,
    #[serde(rename = "monte_carlo")]
    MonteCarlo,
}

#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema, Default)]
#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
pub enum CoordinateSystem {
    #[serde(rename = "cartesian")]
    Cartesian,
    #[serde(rename = "spherical")]
    #[default]
    Spherical,
    #[serde(rename = "hyperspherical")]
    HyperSpherical,
    #[serde(rename = "hyperspherical_flat")]
    HyperSphericalFlat,
    #[serde(rename = "momentum_space")]
    MomentumSpace,
    #[serde(rename = "tropical")]
    MomTrop,
    #[serde(rename = "relative_spherical")]
    RelativeSpherical,
    #[serde(rename = "spherical_common_radial")]
    SphericalCommonRadial,
    #[serde(rename = "spherical_product_common_radial")]
    SphericalProductCommonRadial,
}

impl Default for SamplingSettings {
    fn default() -> Self {
        Self::Default(ParameterizationSettings::default())
    }
}

fn validate_lmb_basis_ids(lmb_basis_ids: &BTreeMap<String, Vec<usize>>) -> Result<(), String> {
    for (graph_name, basis_ids) in lmb_basis_ids {
        if basis_ids.is_empty() {
            return Err(format!(
                "Invalid sampling settings: lmb_basis_ids entry for graph '{graph_name}' cannot be empty."
            ));
        }

        let mut seen = BTreeSet::new();
        for basis_id in basis_ids {
            if !seen.insert(*basis_id) {
                return Err(format!(
                    "Invalid sampling settings: lmb_basis_ids entry for graph '{graph_name}' contains duplicate basis id {basis_id}."
                ));
            }
        }
    }

    Ok(())
}

fn validate_sampling_channel_selection(
    default_selection: &[String],
    selection: &BTreeMap<String, Vec<String>>,
    definitions: &BTreeMap<String, BTreeMap<String, SamplingChannelDefinition>>,
) -> Result<(), String> {
    for (label, entries) in [("default_channel_selection", default_selection)] {
        if entries.iter().any(|entry| entry.trim().is_empty()) {
            return Err(format!(
                "Invalid sampling settings: {label} contains an empty channel selector."
            ));
        }
    }
    for (graph, entries) in selection {
        if graph.trim().is_empty() {
            return Err(
                "Invalid sampling settings: channel_selection contains an empty graph name."
                    .to_string(),
            );
        }
        if entries.iter().any(|entry| entry.trim().is_empty()) {
            return Err(format!(
                "Invalid sampling settings: channel_selection entry for graph '{graph}' contains an empty selector."
            ));
        }
    }
    for (graph, graph_definitions) in definitions {
        if graph.trim().is_empty() {
            return Err(
                "Invalid sampling settings: channel_definitions contains an empty graph name."
                    .to_string(),
            );
        }
        for (name, definition) in graph_definitions {
            if name.trim().is_empty() || definition.around.trim().is_empty() {
                return Err(format!(
                    "Invalid sampling settings: channel definition '{graph}.{name}' requires a non-empty around expression."
                ));
            }
            if definition.parent_lmb.is_empty() {
                return Err(format!(
                    "Invalid sampling settings: channel definition '{graph}.{name}' requires a non-empty parent_lmb."
                ));
            }
            let mut seen = BTreeSet::new();
            if definition.parent_lmb.iter().any(|edge| !seen.insert(edge)) {
                return Err(format!(
                    "Invalid sampling settings: channel definition '{graph}.{name}' contains duplicate parent_lmb edge ids."
                ));
            }
            let mut seen_subspace = BTreeSet::new();
            if definition
                .subspace_lmb
                .iter()
                .any(|edge| !seen_subspace.insert(edge))
            {
                return Err(format!(
                    "Invalid sampling settings: channel definition '{graph}.{name}' contains duplicate subspace_lmb edge ids."
                ));
            }
            let mut seen_on_cut = BTreeSet::new();
            if definition
                .on_cut
                .iter()
                .any(|edge| !seen_on_cut.insert(edge))
            {
                return Err(format!(
                    "Invalid sampling settings: channel definition '{graph}.{name}' contains duplicate on_cut edge ids."
                ));
            }
            SamplingMapDefinition::parse(&definition.around).map_err(|error| {
                format!(
                    "Invalid sampling settings: channel definition '{graph}.{name}' has an invalid around expression: {error}"
                )
            })?;
        }
    }
    Ok(())
}

impl SamplingSettings {
    pub fn as_parser(&self) -> SamplingSettingsParser {
        let mut parser = SamplingSettingsParser::default();
        match self {
            Self::Default(_) => {}
            Self::MultiChanneling(_) => parser.sampling_multichanneling = true,
            Self::DiscreteGraphs(settings) => {
                parser.graphs = SumMode::MonteCarlo;
                parser.graph_names = settings.graph_names.clone();
                parser.orientations = if settings.sample_orientations {
                    SumMode::MonteCarlo
                } else {
                    SumMode::Summed
                };
                parser.sampling_multichanneling = matches!(
                    settings.sampling_type,
                    DiscreteGraphSamplingType::MultiChanneling(_)
                        | DiscreteGraphSamplingType::SamplingMultiChanneling(_)
                );
                if matches!(
                    settings.sampling_type,
                    DiscreteGraphSamplingType::SamplingMultiChanneling(_)
                ) {
                    parser.sampling_channels = SumMode::MonteCarlo;
                }
            }
        }
        match self.get_parameterization_settings() {
            Some(settings) => SamplingSettingsParser {
                alpha: settings.sampling_channels.alpha,
                sampling_channel_weight: settings.sampling_channels.weight,
                coordinate_system: CoordinateSystem::from_mode(settings.mode),
                mapping: settings.mapping,
                b: settings.b,
                power: settings.power,
                lmb_basis_ids: settings.lmb_basis_ids,
                default_channel_selection: settings.sampling_channels.default_channel_selection,
                channel_selection: settings.sampling_channels.channel_selection,
                channel_definitions: settings.sampling_channels.channel_definitions,
                ..parser
            },
            None => SamplingSettingsParser {
                coordinate_system: CoordinateSystem::MomTrop,
                ..parser
            },
        }
    }

    fn from_parser(parser: SamplingSettingsParser) -> Result<Self, String> {
        let SamplingSettingsParser {
            graphs,
            graph_names,
            orientations,
            sampling_multichanneling,
            sampling_channels,
            alpha,
            sampling_channel_weight,
            coordinate_system,
            mapping,
            b,
            power,
            lmb_basis_ids,
            default_channel_selection,
            channel_selection,
            channel_definitions,
        } = parser;

        validate_lmb_basis_ids(&lmb_basis_ids)?;
        let default_channel_selection =
            if sampling_multichanneling && default_channel_selection.is_empty() {
                vec!["auto:optimized_lmb".to_string()]
            } else {
                default_channel_selection
            };
        validate_sampling_channel_selection(
            &default_channel_selection,
            &channel_selection,
            &channel_definitions,
        )?;
        if !alpha.is_finite() || alpha < 0.0 {
            return Err(
                "Invalid sampling settings: alpha must be finite and nonnegative.".to_string(),
            );
        }

        let mut seen_graph_names = BTreeSet::new();
        for graph_name in &graph_names {
            if !seen_graph_names.insert(graph_name) {
                return Err(format!(
                    "Invalid sampling settings: graph_names contains duplicate graph name '{graph_name}'."
                ));
            }
        }
        if !graph_names.is_empty() && matches!(graphs, SumMode::Summed) {
            return Err(
                "Invalid sampling settings: graph_names requires graphs = 'monte_carlo'."
                    .to_string(),
            );
        }

        let sample_orientations = match (graphs.clone(), orientations) {
            (SumMode::Summed, SumMode::Summed) => false,
            (SumMode::Summed, SumMode::MonteCarlo) => {
                return Err(
                    "Invalid sampling settings: orientations can only be set to 'monte_carlo' when graphs is 'monte_carlo'.".to_string(),
                )
            }
            (SumMode::MonteCarlo, SumMode::Summed) => false,
            (SumMode::MonteCarlo, SumMode::MonteCarlo) => true,
        };

        if matches!(coordinate_system, CoordinateSystem::MomTrop) {
            if sampling_channel_weight != SamplingChannelWeight::MapDensity
                || !default_channel_selection.is_empty()
                || !channel_selection.is_empty()
                || !channel_definitions.is_empty()
            {
                return Err(
                    "Invalid sampling settings: advanced sampling channels are incompatible with tropical coordinates."
                        .to_string(),
                );
            }
            if !lmb_basis_ids.is_empty() {
                return Err(
                    "Invalid sampling settings: coordinate_system = 'tropical' is incompatible with lmb_basis_ids."
                        .to_string(),
                );
            }
            if !matches!(graphs, SumMode::MonteCarlo) {
                return Err(
                    "Invalid sampling settings: coordinate_system = 'tropical' requires graphs = 'monte_carlo'."
                        .to_string(),
                );
            }

            if sampling_multichanneling {
                return Err(
                    "Invalid sampling settings: coordinate_system = 'tropical' is incompatible with sampling_multichanneling = true."
                        .to_string(),
                );
            }

            return Ok(SamplingSettings::DiscreteGraphs(
                DiscreteGraphSamplingSettings {
                    graph_names,
                    sample_orientations,
                    sampling_type: DiscreteGraphSamplingType::TropicalSampling(
                        GammaloopTropicalSamplingSettings::default(),
                    ),
                },
            ));
        }

        let mode = coordinate_system.into_mode();
        if matches!(mode, ParameterizationMode::SphericalProductCommonRadial)
            && (!sampling_multichanneling
                || sampling_channel_weight != SamplingChannelWeight::InverseJacobian)
        {
            return Err(
                "Invalid sampling settings: coordinate_system = 'spherical_product_common_radial' requires sampling_multichanneling = true and sampling_channel_weight = 'inverse_jacobian'."
                    .to_string(),
            );
        }
        if sampling_channel_weight == SamplingChannelWeight::InverseJacobian
            && matches!(mode, ParameterizationMode::HyperSphericalFlat)
        {
            return Err(
                "Invalid sampling settings: sampling_channel_weight = 'inverse_jacobian' is incompatible with coordinate_system = 'hyperspherical_flat' because the inverse map is not available."
                    .to_string(),
            );
        }
        if matches!(mapping, ParameterizationMapping::Power) {
            if !matches!(
                mode,
                ParameterizationMode::Spherical
                    | ParameterizationMode::RelativeSpherical
                    | ParameterizationMode::SphericalCommonRadial
                    | ParameterizationMode::SphericalProductCommonRadial
                    | ParameterizationMode::HyperSpherical
                    | ParameterizationMode::HyperSphericalFlat
            ) {
                return Err(
                    "Invalid sampling settings: mapping = 'power' requires a spherical coordinate system."
                        .to_string(),
                );
            }
            if !power.is_finite() || power < 1.0 {
                return Err(
                    "Invalid sampling settings: mapping = 'power' requires power >= 1.".to_string(),
                );
            }
        }
        let parameterization_settings = ParameterizationSettings {
            mode,
            mapping,
            b,
            power,
            lmb_basis_ids,
            sampling_channels: SamplingChannelSelection {
                weight: sampling_channel_weight,
                alpha,
                default_channel_selection,
                channel_selection,
                channel_definitions,
            },
        };

        match graphs {
            SumMode::Summed => {
                if !matches!(sampling_channels, SumMode::Summed) {
                    return Err(
                        "Invalid sampling settings: sampling_channels = 'monte_carlo' requires graphs = 'monte_carlo'."
                            .to_string(),
                    );
                }

                if sample_orientations {
                    return Err(
                        "Invalid sampling settings: orientations can only be Monte Carlo sampled when graphs are Monte Carlo sampled."
                            .to_string(),
                    );
                }

                if sampling_multichanneling {
                    Ok(SamplingSettings::MultiChanneling(MultiChannelingSettings {
                        parameterization_settings,
                    }))
                } else {
                    Ok(SamplingSettings::Default(parameterization_settings))
                }
            }
            SumMode::MonteCarlo => {
                let sampling_type = if sampling_multichanneling {
                    let settings = MultiChannelingSettings {
                        parameterization_settings,
                    };

                    match sampling_channels {
                        SumMode::Summed => DiscreteGraphSamplingType::MultiChanneling(settings),
                        SumMode::MonteCarlo => {
                            DiscreteGraphSamplingType::SamplingMultiChanneling(settings)
                        }
                    }
                } else {
                    DiscreteGraphSamplingType::Default(parameterization_settings)
                };

                Ok(SamplingSettings::DiscreteGraphs(
                    DiscreteGraphSamplingSettings {
                        graph_names,
                        sample_orientations,
                        sampling_type,
                    },
                ))
            }
        }
    }
}

impl Serialize for SamplingSettings {
    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
    where
        S: Serializer,
    {
        self.as_parser().serialize(serializer)
    }
}

impl<'de> Deserialize<'de> for SamplingSettings {
    fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
    where
        D: Deserializer<'de>,
    {
        let parser = SamplingSettingsParser::deserialize(deserializer)?;
        SamplingSettings::from_parser(parser).map_err(serde::de::Error::custom)
    }
}

impl JsonSchema for SamplingSettings {
    fn schema_name() -> std::borrow::Cow<'static, str> {
        "SamplingSettings".into()
    }

    fn json_schema(generator: &mut schemars::SchemaGenerator) -> schemars::Schema {
        SamplingSettingsParser::json_schema(generator)
    }
}

impl CoordinateSystem {
    fn from_mode(mode: ParameterizationMode) -> Self {
        match mode {
            ParameterizationMode::Cartesian => Self::Cartesian,
            ParameterizationMode::Spherical => Self::Spherical,
            ParameterizationMode::RelativeSpherical => Self::RelativeSpherical,
            ParameterizationMode::SphericalCommonRadial => Self::SphericalCommonRadial,
            ParameterizationMode::SphericalProductCommonRadial => {
                Self::SphericalProductCommonRadial
            }
            ParameterizationMode::HyperSpherical => Self::HyperSpherical,
            ParameterizationMode::HyperSphericalFlat => Self::HyperSphericalFlat,
            ParameterizationMode::MomentumSpace => Self::MomentumSpace,
        }
    }

    fn into_mode(self) -> ParameterizationMode {
        match self {
            CoordinateSystem::Cartesian => ParameterizationMode::Cartesian,
            CoordinateSystem::Spherical => ParameterizationMode::Spherical,
            CoordinateSystem::RelativeSpherical => ParameterizationMode::RelativeSpherical,
            CoordinateSystem::SphericalCommonRadial => ParameterizationMode::SphericalCommonRadial,
            CoordinateSystem::SphericalProductCommonRadial => {
                ParameterizationMode::SphericalProductCommonRadial
            }
            CoordinateSystem::HyperSpherical => ParameterizationMode::HyperSpherical,
            CoordinateSystem::HyperSphericalFlat => ParameterizationMode::HyperSphericalFlat,
            CoordinateSystem::MomentumSpace => ParameterizationMode::MomentumSpace,
            CoordinateSystem::MomTrop => {
                unreachable!("tropical coordinate system has no ParameterizationMode equivalent")
            }
        }
    }
}

impl SamplingSettings {
    /// Whether the selected runtime mode maps canonical sampling channels.
    /// Both top-level and discrete graph sampling use the same channel owner;
    /// keeping the predicate central also guards cut-dependent maps in both forms.
    pub(crate) fn uses_sampling_channels(&self) -> bool {
        match self {
            Self::MultiChanneling(_) => true,
            Self::DiscreteGraphs(settings) => matches!(
                &settings.sampling_type,
                DiscreteGraphSamplingType::MultiChanneling(_)
                    | DiscreteGraphSamplingType::SamplingMultiChanneling(_)
            ),
            Self::Default(_) => false,
        }
    }

    pub fn selected_graph_names(&self) -> &[String] {
        match self {
            SamplingSettings::DiscreteGraphs(settings) => &settings.graph_names,
            SamplingSettings::Default(_) | SamplingSettings::MultiChanneling(_) => &[],
        }
    }

    pub fn get_parameterization_settings(&self) -> Option<ParameterizationSettings> {
        match self {
            SamplingSettings::Default(settings) => Some(settings.clone()),
            SamplingSettings::MultiChanneling(settings) => {
                Some(settings.parameterization_settings.clone())
            }
            SamplingSettings::DiscreteGraphs(settings) => match &settings.sampling_type {
                DiscreteGraphSamplingType::Default(settings) => Some(settings.clone()),
                DiscreteGraphSamplingType::MultiChanneling(settings) => {
                    Some(settings.parameterization_settings.clone())
                }
                DiscreteGraphSamplingType::SamplingMultiChanneling(settings) => {
                    Some(settings.parameterization_settings.clone())
                }
                DiscreteGraphSamplingType::TropicalSampling(_) => None,
            },
        }
    }

    pub(crate) fn discrete_depth(&self) -> usize {
        match self {
            SamplingSettings::Default(_) => 0,
            SamplingSettings::MultiChanneling(_) => 0,
            SamplingSettings::DiscreteGraphs(settings) => {
                let depth_from_orientations = settings.sample_orientations as usize;

                match &settings.sampling_type {
                    DiscreteGraphSamplingType::Default(_) => 1 + depth_from_orientations,
                    DiscreteGraphSamplingType::MultiChanneling(_) => 1 + depth_from_orientations,
                    DiscreteGraphSamplingType::SamplingMultiChanneling(_) => {
                        2 + depth_from_orientations
                    }
                    DiscreteGraphSamplingType::TropicalSampling(_) => 1 + depth_from_orientations,
                }
            }
        }
    }

    pub(crate) fn describe_settings(&self) -> String {
        match self {
            SamplingSettings::Default(settings) => {
                format!("{} coordinates", settings.mode)
            }
            SamplingSettings::MultiChanneling(settings) => {
                format!(
                    "lmb multichanneling in {} coordinates",
                    settings.parameterization_settings.mode
                )
            }
            SamplingSettings::DiscreteGraphs(settings) => {
                let discrete_graph_string = if settings.graph_names.is_empty() {
                    "Monte Carlo over graphs".to_string()
                } else {
                    format!(
                        "Monte Carlo over selected graph groups [{}]",
                        settings.graph_names.join(", ")
                    )
                };
                let orientation_sampling_string = if settings.sample_orientations {
                    "and Monte Carlo over orientations"
                } else {
                    ""
                };

                match &settings.sampling_type {
                    DiscreteGraphSamplingType::Default(settings) => {
                        format!(
                            "{} {} in {} coordinates",
                            discrete_graph_string, orientation_sampling_string, settings.mode
                        )
                    }
                    DiscreteGraphSamplingType::MultiChanneling(settings) => {
                        format!(
                            "{}, lmb multichanneling in {} coordinates {}",
                            discrete_graph_string,
                            settings.parameterization_settings.mode,
                            orientation_sampling_string,
                        )
                    }
                    DiscreteGraphSamplingType::SamplingMultiChanneling(settings) => {
                        format!(
                            "{}, {} and monte carlo over lmbs in {} coordinates",
                            discrete_graph_string,
                            orientation_sampling_string,
                            settings.parameterization_settings.mode
                        )
                    }
                    DiscreteGraphSamplingType::TropicalSampling(_) => {
                        format!(
                            "{} {} using 🌴🥥 tropical sampling 🥥🌴",
                            discrete_graph_string, orientation_sampling_string,
                        )
                    }
                }
            }
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Default, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[serde(deny_unknown_fields, default)]
pub struct MultiChannelingSettings {
    /// Continuous loop-momentum parameterization used inside each channel.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub parameterization_settings: ParameterizationSettings,
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[serde(deny_unknown_fields, default)]
pub struct GammaloopTropicalSamplingSettings {
    /// Retry unstable tropical matrix operations at a higher numerical precision.
    #[serde(skip_serializing_if = "is_true")]
    pub upcast_on_failure: bool,
    /// Optional condition threshold used to reject an unstable tropical sampling matrix.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub matrix_stability_test: Option<f64>,
}

impl Default for GammaloopTropicalSamplingSettings {
    fn default() -> Self {
        Self {
            upcast_on_failure: true,
            matrix_stability_test: None,
        }
    }
}

impl GammaloopTropicalSamplingSettings {
    pub fn into_tropical_sampling_settings(&self) -> momtrop::TropicalSamplingSettings {
        momtrop::TropicalSamplingSettings {
            matrix_stability_test: self.matrix_stability_test,
            print_debug_info: false,
            return_metadata: false,
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[serde(tag = "subtype", deny_unknown_fields)]
pub enum DiscreteGraphSamplingType {
    #[serde(rename = "default")]
    Default(ParameterizationSettings),
    #[serde(rename = "multi_channeling")]
    MultiChanneling(MultiChannelingSettings),
    #[serde(rename = "discrete_multi_channeling")]
    SamplingMultiChanneling(MultiChannelingSettings),
    #[serde(rename = "tropical")]
    TropicalSampling(GammaloopTropicalSamplingSettings),
}

impl Default for DiscreteGraphSamplingType {
    fn default() -> Self {
        DiscreteGraphSamplingType::Default(ParameterizationSettings::default())
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema, Default)]
#[serde(deny_unknown_fields, default)]
pub struct DiscreteGraphSamplingSettings {
    /// Optional graph-name allow-list; an empty list includes all generated graphs.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub graph_names: Vec<String>,
    /// Sample graph orientations as an additional discrete integration coordinate.
    #[serde(skip_serializing_if = "is_false")]
    pub sample_orientations: bool,
    /// Continuous and multichannel strategy used inside each selected graph.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub sampling_type: DiscreteGraphSamplingType,
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Default, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[serde(deny_unknown_fields, default)]
pub struct LocalCounterTermSettings {
    /// Widths and override controls for ultraviolet-counterterm localization dampers.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub uv_localisation: UVLocalisationSettings,
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[serde(default, deny_unknown_fields)]
pub struct UVLocalisationSettings {
    /// Width of the transition region used by the ultraviolet sliver damper.
    #[serde(skip_serializing_if = "is_float::<10>")]
    pub sliver_width: f64,
    /// Make the compact sliver profile smooth at its support boundary.
    #[serde(skip_serializing_if = "is_false")]
    pub smooth_sliver: bool,
    /// Derive the localization width dynamically from the sampled ultraviolet configuration.
    #[serde(skip_serializing_if = "is_false")]
    pub dynamic_width: bool,
    /// Width of the Gaussian ultraviolet localization factor.
    #[serde(skip_serializing_if = "is_float::<1>")]
    pub gaussian_width: f64,
    /// Replace ultraviolet localization dampers by one for diagnostic comparisons.
    #[serde(skip_serializing_if = "is_false")]
    pub force_uv_dampers_to_one: bool,
}

impl Default for UVLocalisationSettings {
    fn default() -> Self {
        Self {
            sliver_width: 10.0,
            smooth_sliver: false,
            dynamic_width: false,
            gaussian_width: 1.0,
            force_uv_dampers_to_one: false,
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Default, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[serde(deny_unknown_fields, default)]
pub struct IntegratedCounterTermSettings {
    /// Integration domain and damping function used for integrated counterterms.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub range: IntegratedCounterTermRange,
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[serde(tag = "type")]
#[serde(deny_unknown_fields)]
pub enum IntegratedCounterTermRange {
    /// Integrate over the full radial domain and regulate it with an H-function profile.
    #[serde(rename = "infinite")]
    Infinite {
        /// Damping profile used to make the infinite-domain counterterm integral convergent.
        h_function_settings: HFunctionSettings,
    },
    /// Integrate over the compact radial domain without an infinite-tail damping profile.
    #[serde(rename = "compact")]
    Compact {},
}

impl Default for IntegratedCounterTermRange {
    fn default() -> Self {
        Self::Infinite {
            h_function_settings: HFunctionSettings::default(),
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object, get_all, set_all))]
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[serde(deny_unknown_fields)]
#[serde(default)]
pub struct OverlapSettings {
    //v
    /// Explicit three-vector center used for every threshold overlap, when provided.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub force_global_center: Option<Vec<[f64; 3]>>,
    /// Validate that the selected global center lies inside every required E-surface.
    #[serde(skip_serializing_if = "is_true")]
    pub check_global_center: bool,
    /// Try the loop-momentum origin as the global center before solving for another point.
    #[serde(skip_serializing_if = "is_true")]
    pub try_origin: bool,
    /// Try the origin independently for every available loop-momentum basis.
    #[serde(skip_serializing_if = "is_false")]
    pub try_origin_all_lmbs: bool,
}

impl Default for OverlapSettings {
    fn default() -> Self {
        Self {
            force_global_center: None,
            check_global_center: true,
            try_origin: true,
            try_origin_all_lmbs: false,
        }
    }
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object))]
#[derive(Debug, Clone, Default, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[serde(deny_unknown_fields)]
// #[trait_decode(trait= GammaLoopContext)]
pub enum HFunction {
    #[default]
    #[serde(rename = "poly_exponential")]
    PolyExponential,
    #[serde(rename = "exponential")]
    Exponential,
    #[serde(rename = "poly_left_right_exponential")]
    PolyLeftRightExponential,
    #[serde(rename = "exponential_ct")]
    ExponentialCT,
}

#[cfg_attr(feature = "python", pyo3::pyclass(from_py_object))]
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[serde(deny_unknown_fields)]
#[serde(default)]
pub struct HFunctionSettings {
    /// Functional form of the counterterm damping profile.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub function: HFunction,
    /// Characteristic width of the damping profile.
    #[serde(skip_serializing_if = "is_float::<1>")]
    pub sigma: f64,
    /// Apply the damping profile; disabling it leaves the corresponding factor equal to one.
    #[serde(skip_serializing_if = "is_true")]
    pub enabled_dampening: bool,
    /// Optional polynomial power overriding the profile's built-in choice.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub power: Option<usize>,
}

impl Default for HFunctionSettings {
    fn default() -> Self {
        Self {
            sigma: 1.0,
            function: HFunction::default(),
            enabled_dampening: true,
            power: None,
        }
    }
}

impl HFunctionSettings {
    /// Supported powers of the polynomial profiles.
    const NORMALIZED_POWERS: [usize; 12] = [0, 1, 3, 4, 6, 7, 9, 10, 12, 13, 15, 16];

    /// Integral of the unnormalized profile over `t > 0`, in the precision of
    /// `sigma`. With `x = t/sigma`, `exponential` integrates `exp(-x^2)` to
    /// `sqrt(pi)/2`, while the polynomial profiles are Bessel integrals:
    /// `x^-p exp(2 - x^2 - x^-2)` gives `e^2 K_{(1-p)/2}(2)` and
    /// `x^-p exp(2 - x - 1/x)` gives `2 e^2 K_{1-p}(2)`.
    pub(crate) fn normalization<T: FloatLike>(&self, sigma: &F<T>) -> F<T> {
        let one = sigma.one();
        let power = self.power.unwrap_or(0);
        if self.function != HFunction::Exponential && !Self::NORMALIZED_POWERS.contains(&power) {
            panic!("Value {power} of power in poly exponential h function not supported");
        }
        let unit = match self.function {
            HFunction::Exponential => Self::scaled_bessel_k_at_two(1, &one),
            HFunction::PolyExponential => Self::scaled_bessel_k_at_two(power.abs_diff(1), &one),
            HFunction::PolyLeftRightExponential => {
                one.from_usize(2) * Self::scaled_bessel_k_at_two(2 * power.abs_diff(1), &one)
            }
            HFunction::ExponentialCT => {
                unreachable!("exponential_ct is a threshold localization without a normalization")
            }
        };
        unit * sigma
    }

    /// `e^2 K_nu(2)` for `nu = twice_order/2`, in the precision of `one`.
    ///
    /// Half-integer orders are elementary, `e^2 K_{n+1/2}(2) = sqrt(pi)/2 *
    /// sum_k (n+k)!/(k! (n-k)! 4^k)`, an exact rational times `sqrt(pi)/2`.
    /// Integer orders start from the `K_0` series, where `log(z/2)` vanishes at
    /// `z = 2` and leaves Euler's constant as the only transcendental input.
    /// The Wronskian `I_0 K_1 + I_1 K_0 = 1/2` gives `K_1`, and the upward
    /// recurrence `K_{n+1}(2) = K_{n-1}(2) + n K_n(2)` is stable for `K`.
    fn scaled_bessel_k_at_two<T: FloatLike>(twice_order: usize, one: &F<T>) -> F<T> {
        use std::{
            any::{Any, TypeId},
            cell::RefCell,
            collections::HashMap,
        };
        use symbolica::domains::float::SingleFloat;

        if twice_order % 2 == 1 {
            let n = twice_order / 2;
            let sum = (0..=n).fold(Rational::zero(), |sum, k| {
                sum + Rational::from((
                    Integer::factorial((n + k) as u32),
                    Integer::factorial(k as u32)
                        * Integer::factorial((n - k) as u32)
                        * Integer::from(4).pow(k as u64),
                ))
            });
            return one.PI().sqrt() * one.from_rational(&(sum / Rational::from(2)));
        }
        // Cache native unit values per worker, like the decimal-promotion caches.
        // Include both the numeric representation and its precision; sigma is
        // applied by normalization() and must never become part of this value.
        type NativeValues = HashMap<(TypeId, u32, usize), Box<dyn Any>>;
        thread_local! {
            static CACHE: RefCell<NativeValues> = RefCell::default();
        }
        let key = (TypeId::of::<T>(), one.0.get_precision(), twice_order);
        if let Some(value) = CACHE.with_borrow(|cache| {
            cache.get(&key).map(|value| {
                value
                    .downcast_ref::<F<T>>()
                    .expect("cache key includes the numeric type")
                    .clone()
            })
        }) {
            return value;
        }

        // Each series term carries (z^2/4)^k/(k!)^2 = 1/(k!)^2 at z = 2.
        let (mut i0, mut i1, mut harmonic_sum) = (one.zero(), one.zero(), one.zero());
        let (mut inverse_square, mut harmonic) = (one.clone(), one.zero());
        for k in 1.. {
            i0 += &inverse_square;
            i1 += &inverse_square / one.from_usize(k);
            inverse_square = &inverse_square / one.from_usize(k * k);
            harmonic += one.from_usize(k).inv();
            // H_k >= 1 bounds every remaining I_0 and I_1 term as well.
            let term = &harmonic * &inverse_square;
            harmonic_sum += &term;
            if term <= &harmonic_sum * one.epsilon() {
                break;
            }
        }
        let k0 = harmonic_sum - one.euler() * &i0;
        let k1 = (one.from_usize(2).inv() - &i1 * &k0) / &i0;
        let order = twice_order / 2;
        let (mut lower, mut upper) = (k0, k1);
        for n in 1..order {
            (lower, upper) = (upper.clone(), lower + one.from_usize(n) * &upper);
        }
        let value = one.from_usize(2).exp() * if order == 0 { lower } else { upper };
        CACHE.with_borrow_mut(|cache| {
            // Bound retention for callers using many numeric types or precisions.
            if cache.len() == 1024 {
                cache.clear();
            }
            cache.insert(key, Box::new(value.clone()));
        });
        value
    }

    /// Require a profile normalized to one over positive scales `t`. LU-h
    /// sampling and Fermi-surface localization both insert it as a unit integral.
    pub(crate) fn validate_normalized_profile(&self) -> Result<()> {
        if !self.sigma.is_finite() || self.sigma <= 0.0 {
            return Err(eyre::eyre!(
                "A normalized h_function profile requires a positive finite h_function.sigma; got {}",
                self.sigma
            ));
        }
        match self.function {
            HFunction::Exponential => Ok(()),
            HFunction::PolyExponential | HFunction::PolyLeftRightExponential => {
                let power = self.power.unwrap_or(0);
                if Self::NORMALIZED_POWERS.contains(&power) {
                    Ok(())
                } else {
                    Err(eyre::eyre!(
                        "Unsupported h_function.power {power}; supported powers are {:?}",
                        Self::NORMALIZED_POWERS
                    ))
                }
            }
            HFunction::ExponentialCT => Err(eyre::eyre!(
                "exponential_ct is a local threshold localization, not a normalized h_function profile"
            )),
        }
    }
}