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::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;

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

#[cfg_attr(
    feature = "python_api",
    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_api",
    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_api", 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_api",
    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,
    /// 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,

            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_api", 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_api",
    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_api",
    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_api",
    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>>,
}

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(),
        }
    }
}

#[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 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,
}

#[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 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 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_api",
    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 {
    /// 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_true")]
    pub escalate_if_exact_zero: bool,
    /// 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 {
            rotation_axis: _default_rotation_axis(),
            levels: _default_stability_levels(),
            check_on_norm: true,
            escalate_if_exact_zero: false,
            loop_momenta_norm_escalation_factor: -1.0,
            recording: None,
        }
    }
}

#[cfg_attr(
    feature = "python_api",
    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_api",
    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,
    /// 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,
            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,
            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,
            escalate_for_large_weight_threshold: -1.0,
        }
    }
}

#[cfg_attr(
    feature = "python_api",
    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_api", 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_api",
    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_api",
    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_api", 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_api",
    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_api",
    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 loop-momentum bases.
    #[serde(skip_serializing_if = "is_false")]
    pub lmb_multichanneling: bool,
    /// Whether loop-momentum-basis channels are summed or sampled discretely.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub lmb_channels: SumMode,
    /// Exponent controlling the sharpness of multichannel weights.
    #[serde(skip_serializing_if = "is_float::<3>")]
    pub alpha: f64,
    /// Rule used to construct the relative probability of each loop-momentum-basis channel.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub lmb_channel_weight: LmbChannelWeight,
    /// 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>>,
}

impl Default for SamplingSettingsParser {
    fn default() -> Self {
        Self {
            graphs: SumMode::Summed,
            graph_names: Vec::new(),
            orientations: SumMode::Summed,
            lmb_multichanneling: false,
            lmb_channels: SumMode::Summed,
            alpha: 3.0,
            lmb_channel_weight: LmbChannelWeight::default(),
            coordinate_system: CoordinateSystem::Spherical,
            mapping: ParameterizationMapping::Linear,
            b: 1.0,
            power: 1.0,
            lmb_basis_ids: BTreeMap::new(),
        }
    }
}

#[derive(
    Debug, Default, Clone, Copy, Serialize, Deserialize, PartialEq, Eq, Encode, Decode, JsonSchema,
)]
#[cfg_attr(
    feature = "python_api",
    pyo3::pyclass(from_py_object, get_all, set_all)
)]
pub enum LmbChannelWeight {
    #[serde(rename = "ose")]
    #[default]
    Ose,
    #[serde(rename = "inverse_jacobian")]
    InverseJacobian,
}

#[derive(Debug, Default, Clone, Serialize, Deserialize, PartialEq, Encode, Decode, JsonSchema)]
#[cfg_attr(
    feature = "python_api",
    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_api",
    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(())
}

impl SamplingSettings {
    fn as_parser(&self) -> SamplingSettingsParser {
        match self {
            SamplingSettings::Default(settings) => SamplingSettingsParser {
                graphs: SumMode::Summed,
                graph_names: Vec::new(),
                orientations: SumMode::Summed,
                lmb_multichanneling: false,
                lmb_channels: SumMode::Summed,
                alpha: 3.0,
                lmb_channel_weight: LmbChannelWeight::default(),
                coordinate_system: CoordinateSystem::from_mode(settings.mode.clone()),
                mapping: settings.mapping.clone(),
                b: settings.b,
                power: settings.power,
                lmb_basis_ids: settings.lmb_basis_ids.clone(),
            },
            SamplingSettings::MultiChanneling(settings) => SamplingSettingsParser {
                graphs: SumMode::Summed,
                graph_names: Vec::new(),
                orientations: SumMode::Summed,
                lmb_multichanneling: true,
                lmb_channels: SumMode::Summed,
                alpha: settings.alpha,
                lmb_channel_weight: settings.channel_weight,
                coordinate_system: CoordinateSystem::from_mode(
                    settings.parameterization_settings.mode.clone(),
                ),
                mapping: settings.parameterization_settings.mapping.clone(),
                b: settings.parameterization_settings.b,
                power: settings.parameterization_settings.power,
                lmb_basis_ids: settings.parameterization_settings.lmb_basis_ids.clone(),
            },
            SamplingSettings::DiscreteGraphs(settings) => {
                let orientations = if settings.sample_orientations {
                    SumMode::MonteCarlo
                } else {
                    SumMode::Summed
                };

                match &settings.sampling_type {
                    DiscreteGraphSamplingType::Default(parameterization_settings) => {
                        SamplingSettingsParser {
                            graphs: SumMode::MonteCarlo,
                            graph_names: settings.graph_names.clone(),
                            orientations,
                            lmb_multichanneling: false,
                            lmb_channels: SumMode::Summed,
                            alpha: 3.0,
                            lmb_channel_weight: LmbChannelWeight::default(),
                            coordinate_system: CoordinateSystem::from_mode(
                                parameterization_settings.mode.clone(),
                            ),
                            mapping: parameterization_settings.mapping.clone(),
                            b: parameterization_settings.b,
                            power: parameterization_settings.power,
                            lmb_basis_ids: parameterization_settings.lmb_basis_ids.clone(),
                        }
                    }
                    DiscreteGraphSamplingType::MultiChanneling(multichanneling_settings) => {
                        SamplingSettingsParser {
                            graphs: SumMode::MonteCarlo,
                            graph_names: settings.graph_names.clone(),
                            orientations,
                            lmb_multichanneling: true,
                            lmb_channels: SumMode::Summed,
                            alpha: multichanneling_settings.alpha,
                            lmb_channel_weight: multichanneling_settings.channel_weight,
                            coordinate_system: CoordinateSystem::from_mode(
                                multichanneling_settings
                                    .parameterization_settings
                                    .mode
                                    .clone(),
                            ),
                            mapping: multichanneling_settings
                                .parameterization_settings
                                .mapping
                                .clone(),
                            b: multichanneling_settings.parameterization_settings.b,
                            power: multichanneling_settings.parameterization_settings.power,
                            lmb_basis_ids: multichanneling_settings
                                .parameterization_settings
                                .lmb_basis_ids
                                .clone(),
                        }
                    }
                    DiscreteGraphSamplingType::DiscreteMultiChanneling(
                        multichanneling_settings,
                    ) => SamplingSettingsParser {
                        graphs: SumMode::MonteCarlo,
                        graph_names: settings.graph_names.clone(),
                        orientations,
                        lmb_multichanneling: true,
                        lmb_channels: SumMode::MonteCarlo,
                        alpha: multichanneling_settings.alpha,
                        lmb_channel_weight: multichanneling_settings.channel_weight,
                        coordinate_system: CoordinateSystem::from_mode(
                            multichanneling_settings
                                .parameterization_settings
                                .mode
                                .clone(),
                        ),
                        mapping: multichanneling_settings
                            .parameterization_settings
                            .mapping
                            .clone(),
                        b: multichanneling_settings.parameterization_settings.b,
                        power: multichanneling_settings.parameterization_settings.power,
                        lmb_basis_ids: multichanneling_settings
                            .parameterization_settings
                            .lmb_basis_ids
                            .clone(),
                    },
                    DiscreteGraphSamplingType::TropicalSampling(_) => SamplingSettingsParser {
                        graphs: SumMode::MonteCarlo,
                        graph_names: settings.graph_names.clone(),
                        orientations,
                        lmb_multichanneling: false,
                        lmb_channels: SumMode::Summed,
                        alpha: 3.0,
                        lmb_channel_weight: LmbChannelWeight::default(),
                        coordinate_system: CoordinateSystem::MomTrop,
                        mapping: ParameterizationMapping::default(),
                        b: 1.0,
                        power: 1.0,
                        lmb_basis_ids: BTreeMap::new(),
                    },
                }
            }
        }
    }

    fn from_parser(parser: SamplingSettingsParser) -> Result<Self, String> {
        let SamplingSettingsParser {
            graphs,
            graph_names,
            orientations,
            lmb_multichanneling,
            lmb_channels,
            alpha,
            lmb_channel_weight,
            coordinate_system,
            mapping,
            b,
            power,
            lmb_basis_ids,
        } = parser;

        validate_lmb_basis_ids(&lmb_basis_ids)?;

        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 !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 lmb_multichanneling {
                return Err(
                    "Invalid sampling settings: coordinate_system = 'tropical' is incompatible with lmb_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)
            && (!lmb_multichanneling || lmb_channel_weight != LmbChannelWeight::InverseJacobian)
        {
            return Err(
                "Invalid sampling settings: coordinate_system = 'spherical_product_common_radial' requires lmb_multichanneling = true and lmb_channel_weight = 'inverse_jacobian'."
                    .to_string(),
            );
        }
        if lmb_channel_weight == LmbChannelWeight::InverseJacobian
            && matches!(mode, ParameterizationMode::HyperSphericalFlat)
        {
            return Err(
                "Invalid sampling settings: lmb_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,
        };

        match graphs {
            SumMode::Summed => {
                if !matches!(lmb_channels, SumMode::Summed) {
                    return Err(
                        "Invalid sampling settings: lmb_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 lmb_multichanneling {
                    Ok(SamplingSettings::MultiChanneling(MultiChannelingSettings {
                        alpha,
                        channel_weight: lmb_channel_weight,
                        parameterization_settings,
                    }))
                } else {
                    Ok(SamplingSettings::Default(parameterization_settings))
                }
            }
            SumMode::MonteCarlo => {
                let sampling_type = if lmb_multichanneling {
                    let settings = MultiChannelingSettings {
                        alpha,
                        channel_weight: lmb_channel_weight,
                        parameterization_settings,
                    };

                    match lmb_channels {
                        SumMode::Summed => DiscreteGraphSamplingType::MultiChanneling(settings),
                        SumMode::MonteCarlo => {
                            DiscreteGraphSamplingType::DiscreteMultiChanneling(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 {
    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::DiscreteMultiChanneling(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::DiscreteMultiChanneling(_) => {
                        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::DiscreteMultiChanneling(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_api",
    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 MultiChannelingSettings {
    /// Exponent controlling the sharpness of loop-momentum-basis channel weights.
    #[serde(skip_serializing_if = "is_float::<3>")]
    pub alpha: f64,
    /// Rule used to assign relative loop-momentum-basis channel weights.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub channel_weight: LmbChannelWeight,
    /// Continuous loop-momentum parameterization used inside each channel.
    #[serde(skip_serializing_if = "IsDefault::is_default")]
    pub parameterization_settings: ParameterizationSettings,
}

impl Default for MultiChannelingSettings {
    fn default() -> Self {
        Self {
            alpha: 3.0,
            channel_weight: LmbChannelWeight::default(),
            parameterization_settings: ParameterizationSettings::default(),
        }
    }
}

#[cfg_attr(
    feature = "python_api",
    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_api",
    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")]
    DiscreteMultiChanneling(MultiChannelingSettings),
    #[serde(rename = "tropical")]
    TropicalSampling(GammaloopTropicalSamplingSettings),
}

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

#[cfg_attr(
    feature = "python_api",
    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_api",
    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_api",
    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,
    /// 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,
            dynamic_width: false,
            gaussian_width: 1.0,
            force_uv_dampers_to_one: false,
        }
    }
}

#[cfg_attr(
    feature = "python_api",
    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_api",
    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_api",
    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_api", 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_api", 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,
        }
    }
}