gammalooprs/processes/

cross_section.rs

use std::{
    collections::{BTreeMap, BTreeSet},
    fmt::{Display, Formatter},
    fs::{self, File},
    io::Write,
    ops::{Index, IndexMut},
    path::{Path, PathBuf},
    sync::Mutex,
};

#[cfg(test)]
mod gauge_sewing_tests;
#[cfg(test)]
mod sewing_tests;

// use bincode::{Decode, Encode};
use bincode_trait_derive::{Decode, Encode};
use color_eyre::Result;
use itertools::Itertools;
use rayon::{
    ThreadPool,
    iter::{IntoParallelRefMutIterator, ParallelIterator},
};
use spenso::algebra::algebraic_traits::IsZero;
use tracing::info;
use vakint::Vakint;

use crate::{
    DependentMomentaConstructor, GammaLoopContext, GammaLoopContextContainer,
    cff::{
        CffEnergyDegreeBoundReport, CutCFFIndex,
        esurface::{RaisedEsurfaceData, RaisedEsurfaceGroup, RaisedEsurfaceId},
        expression::{
            OrientationID, normalize_cut_edge_support_with_raised_edge_groups,
            normalize_three_d_expression_cut_support_with_raised_edge_groups,
        },
        surface::HybridSurfaceID,
    },
    debug_tags, define_index,
    feyngen::GenerationType,
    graph::{
        GraphGroup, GroupId, LMBext, LmbChannelFallback, LmbIndex, LoopMomentumBasis,
        ThresholdPinchStatus,
        cuts::{CutSet, LuCutSelection, ResidueSelector},
        edge::EdgeMass,
        parse::complete_group_parsing,
        threshold_counterterms::{ThresholdCountertermMultiplier, ThresholdCountertermVariant},
    },
    integrands::process::{
        GenericEvaluator, LmbMultiChannelingSetup, ParamBuilder,
        cross_section::CrossSectionIntegrandData, graph_to_group_id_for_group_structure,
    },
    model::ArcParticle,
    momentum::{
        Helicity,
        sample::{ExternalIndex, SubspaceData},
    },
    processes::{
        DotExportSettings, EvaluatorSettings, GraphCutSelectionSubject, GraphGenerationStats,
        GraphGroupSelectionMode, GraphGroupSelectionPlan, GraphGroupSelectionSpec,
        GraphSelectionSignatureInventory, GraphSelectionSubject, NamedGraphGenerationReport,
        RaisedCutSignatureInventory,
    },
    settings::{
        GlobalSettings, RuntimeSettings, global::GenerationSettings, runtime::LockedRuntimeSettings,
    },
    utils::{
        DEFAULT_ESURFACE_EXISTENCE_THRESHOLD, F, GS, W_,
        hyperdual_utils::{shape_from_cut_cff_index, simple_n_deriv_shape},
    },
    uv::{
        approx::{CutStructure, OrientationProjection},
        forest::ParametricIntegrands,
    },
};
use eyre::{Context, eyre};
use linnet::half_edge::{
    involution::{EdgeIndex, EdgeVec, Orientation},
    subgraph::{
        HedgeNode, Inclusion, InternalSubGraph, ModifySubSet, OrientedCut, SuBitGraph,
        SubGraphLike, SubSetLike, SubSetOps,
    },
};
use serde::{Deserialize, Serialize};
use symbolica::{domains::dual::HyperDual, prelude::*};
use three_dimensional_reps::GeneratedThreeDExpression;
use tracing::{debug, warn};
use typed_index_collections::{TiVec, ti_vec};

use super::generation_progress::{self, GenerationProcessKind, GenerationProgressPhase};
use super::threshold_counterterms::{
    ResolvedCutGroupThresholdCounterterms, ResolvedThresholdCountertermAssociation,
    ResolvedThresholdCountertermVariant, ResolvedThresholdCounterterms, ThresholdCountertermOrigin,
    ThresholdCountertermSide, ThresholdCountertermVariantId,
};

use crate::{
    cff::esurface::{Esurface, EsurfaceID},
    graph::{ExternalConnection, FeynmanGraph, Graph},
    integrands::process::{
        ProcessIntegrand,
        cross_section::{CrossSectionGraphTerm, CrossSectionIntegrand},
    },
    model::Model,
};

use crate::processes::ProcessDefinition;

#[derive(Clone, Debug, Encode, Decode)]
pub struct IteratedCtCollection<T> {
    data: Vec<T>,
    num_right_thresholds: usize,
}

#[cfg(test)]
#[path = "cross_section/ir_safe_threshold_fixture_tests.rs"]
mod ir_safe_threshold_fixture_tests;

impl<T> IteratedCtCollection<T> {
    pub(crate) fn new(
        data: Vec<T>,
        num_left_thresholds: usize,
        num_right_thresholds: usize,
    ) -> Self {
        let expected_len = num_left_thresholds
            .checked_mul(num_right_thresholds)
            .expect("iterated threshold-counterterm dimensions overflow usize");
        assert_eq!(
            data.len(),
            expected_len,
            "iterated threshold-counterterm data must contain one entry per left/right pair"
        );
        Self {
            data,
            num_right_thresholds,
        }
    }

    pub fn map_ref<U, F>(&self, f: F) -> IteratedCtCollection<U>
    where
        F: Fn(&T) -> U,
    {
        let data = self.data.iter().map(f).collect();
        IteratedCtCollection {
            data,
            num_right_thresholds: self.num_right_thresholds,
        }
    }

    pub(crate) fn iter(&self) -> impl Iterator<Item = &T> {
        self.data.iter()
    }

    pub(crate) fn iter_mut(&mut self) -> impl Iterator<Item = &mut T> {
        self.data.iter_mut()
    }

    pub(crate) fn num_right_thresholds(&self) -> usize {
        self.num_right_thresholds
    }
}

impl<T> Index<(LeftThresholdId, RightThresholdId)> for IteratedCtCollection<T> {
    type Output = T;

    fn index(&self, index: (LeftThresholdId, RightThresholdId)) -> &Self::Output {
        let (left_id, right_id) = index;
        &self.data[left_id.0 * self.num_right_thresholds + right_id.0]
    }
}

impl<T> IndexMut<(LeftThresholdId, RightThresholdId)> for IteratedCtCollection<T> {
    fn index_mut(&mut self, index: (LeftThresholdId, RightThresholdId)) -> &mut Self::Output {
        let (left_id, right_id) = index;
        &mut self.data[left_id.0 * self.num_right_thresholds + right_id.0]
    }
}

#[derive(Clone, Encode, Decode)]
#[trait_decode(trait = GammaLoopContext)]
pub struct LUCounterTermData {
    pub left_thresholds: TiVec<LeftThresholdId, RaisedEsurfaceGroup>,
    pub right_thresholds: TiVec<RightThresholdId, RaisedEsurfaceGroup>,
    pub left_variant_ids: TiVec<LeftThresholdId, ThresholdCountertermVariantId>,
    pub right_variant_ids: TiVec<RightThresholdId, ThresholdCountertermVariantId>,
    pub left_subspaces: TiVec<LeftThresholdId, SubspaceData>,
    pub right_subspaces: TiVec<RightThresholdId, SubspaceData>,
    pub left_atoms: TiVec<LeftThresholdId, ParametricIntegrands>,
    pub right_atoms: TiVec<RightThresholdId, ParametricIntegrands>,
    pub iterated: IteratedCtCollection<ParametricIntegrands>,
}

#[derive(Clone, Debug)]
pub(crate) struct SingleThresholdPieces<T> {
    pub(crate) local: T,
    pub(crate) integrated: T,
}

#[derive(Clone, Debug)]
pub(crate) struct IteratedThresholdPieces<T> {
    pub(crate) local_local: T,
    pub(crate) local_integrated: T,
    pub(crate) integrated_local: T,
    pub(crate) integrated_integrated: T,
}

#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub(crate) enum LUThresholdHelperOutputs {
    Legacy,
    Pieces,
    LegacyAndPieces,
}

fn max_dual_size_for_cut_cff_indices<'a>(
    cut_cff_indices: impl Iterator<Item = &'a CutCFFIndex>,
) -> usize {
    cut_cff_indices
        .map(|cut_cff_index| {
            shape_from_cut_cff_index(cut_cff_index)
                .map(|shape| HyperDual::<F<f64>>::new(shape).values.len())
                .unwrap_or(1)
        })
        .max()
        .unwrap_or(1)
}

define_index! {pub struct RightThresholdId;}
define_index! {pub struct LeftThresholdId;}
define_index! {pub struct CutGroupId;}
define_index! {pub struct TopologicalThresholdId;}

/// Eligibility of one threshold E-surface relative to one side of one physical cut.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum ThresholdCountertermStatus {
    /// The E-surface does not depend on any active loop coordinate of the cut-side subspace.
    NoRadialDependence,
    /// The zero can only occur at a pinched boundary and has no strictly existing region.
    AlwaysPinched,
    /// A strictly existing region is possible, but its boundary includes a pinched limit.
    CanBecomePinched,
    /// Current resolved masses and kinematics rigorously exclude a zero over the full cut phase space.
    ProvenNonExisting,
    /// A zero is possible, or the available metadata is insufficient to exclude one safely.
    PotentiallyExisting,
}

impl ThresholdCountertermStatus {
    pub fn is_eligible_for_generation(self, check_current_model: bool) -> bool {
        match self {
            Self::NoRadialDependence | Self::AlwaysPinched => false,
            Self::ProvenNonExisting => !check_current_model,
            Self::CanBecomePinched | Self::PotentiallyExisting => true,
        }
    }
}

/// Persisted topology and invariant-bound metadata for one generated cut-relative CT association.
#[derive(Clone, Debug, Encode, Decode)]
#[trait_decode(trait = GammaLoopContext)]
pub struct ThresholdCountertermAssociation {
    pub topological_threshold_id: TopologicalThresholdId,
    pub esurface_id: EsurfaceID,
    pub cut_boundary_edges: Vec<EdgeIndex>,
    pub threshold_boundary_edges: Vec<EdgeIndex>,
    pub invariant_bound_is_applicable: bool,
}

#[derive(Clone, Debug, Default, Encode, Decode)]
#[trait_decode(trait = GammaLoopContext)]
pub struct CutThresholdCountertermAssociations {
    pub left: Vec<ThresholdCountertermAssociation>,
    pub right: Vec<ThresholdCountertermAssociation>,
}

impl ThresholdCountertermAssociation {
    fn mass_sum(edges: &[EdgeIndex], masses: &EdgeVec<Option<F<f64>>>) -> Option<F<f64>> {
        edges.iter().try_fold(F(0.0), |sum, edge_id| {
            masses[*edge_id].as_ref().map(|mass| sum + mass.abs())
        })
    }

    fn incoming_invariant_mass(
        graph: &Graph,
        runtime_settings: &RuntimeSettings,
    ) -> Option<F<f64>> {
        let external_momenta = runtime_settings
            .kinematics
            .externals
            .get_dependent_externals(DependentMomentaConstructor::CrossSection)
            .ok()?;
        if external_momenta.is_empty()
            || external_momenta.len() != graph.get_edges_in_initial_state_cut().len()
        {
            return None;
        }

        let total = external_momenta
            .into_iter()
            .reduce(|sum, momentum| sum + momentum)?;
        let invariant_squared = total.square();
        if invariant_squared.is_nan()
            || invariant_squared.is_infinite()
            || invariant_squared <= F(0.0)
        {
            None
        } else {
            Some(invariant_squared.sqrt())
        }
    }

    fn classify_from_invariant_bounds(
        &self,
        cut_boundary_mass: &F<f64>,
        threshold_boundary_mass: &F<f64>,
        maximum_cut_invariant: Option<F<f64>>,
        tolerance: &F<f64>,
    ) -> ThresholdCountertermStatus {
        if !self.invariant_bound_is_applicable
            || cut_boundary_mass.is_nan()
            || cut_boundary_mass.is_infinite()
            || threshold_boundary_mass.is_nan()
            || threshold_boundary_mass.is_infinite()
            || tolerance.is_nan()
            || tolerance.is_infinite()
            || maximum_cut_invariant
                .as_ref()
                .is_some_and(|bound| bound.is_nan() || bound.is_infinite())
        {
            return ThresholdCountertermStatus::PotentiallyExisting;
        }

        let minimum_margin = cut_boundary_mass * cut_boundary_mass
            - threshold_boundary_mass * threshold_boundary_mass;
        let masses_match = minimum_margin.abs() <= *tolerance;
        let threshold_minimum_is_reachable = minimum_margin <= *tolerance;

        if masses_match
            && (self.cut_boundary_edges.len() <= 1 || self.threshold_boundary_edges.len() <= 1)
        {
            return ThresholdCountertermStatus::AlwaysPinched;
        }

        if self.cut_boundary_edges.len() == 1 && self.threshold_boundary_edges.len() == 1 {
            return ThresholdCountertermStatus::ProvenNonExisting;
        }

        if self.cut_boundary_edges.len() > 1 && self.threshold_boundary_edges.len() == 1 {
            let lower_margin = threshold_boundary_mass * threshold_boundary_mass
                - cut_boundary_mass * cut_boundary_mass;
            if lower_margin <= *tolerance {
                return ThresholdCountertermStatus::ProvenNonExisting;
            }
        }

        let Some(maximum_cut_invariant) = maximum_cut_invariant else {
            return if masses_match
                && self.cut_boundary_edges.len() > 1
                && self.threshold_boundary_edges.len() > 1
            {
                ThresholdCountertermStatus::CanBecomePinched
            } else {
                ThresholdCountertermStatus::PotentiallyExisting
            };
        };

        let maximum_margin = maximum_cut_invariant * maximum_cut_invariant
            - threshold_boundary_mass * threshold_boundary_mass;
        if maximum_margin.abs() <= *tolerance {
            ThresholdCountertermStatus::AlwaysPinched
        } else if maximum_margin < -*tolerance {
            ThresholdCountertermStatus::ProvenNonExisting
        } else if threshold_minimum_is_reachable
            && self.cut_boundary_edges.len() > 1
            && self.threshold_boundary_edges.len() > 1
        {
            ThresholdCountertermStatus::CanBecomePinched
        } else {
            ThresholdCountertermStatus::PotentiallyExisting
        }
    }

    pub fn classify_for_model(
        &self,
        graph: &Graph,
        cut: &CrossSectionCut,
        model: &Model,
        param_builder: &ParamBuilder,
        runtime_settings: &RuntimeSettings,
        normalized_margin_tolerance: f64,
    ) -> ThresholdCountertermStatus {
        if !self.invariant_bound_is_applicable {
            return ThresholdCountertermStatus::PotentiallyExisting;
        }

        match graph
            .classify_threshold_pinch(&self.cut_boundary_edges, &self.threshold_boundary_edges)
        {
            ThresholdPinchStatus::Always => {
                return ThresholdCountertermStatus::AlwaysPinched;
            }
            ThresholdPinchStatus::CanBecome | ThresholdPinchStatus::NotProven => {}
        }

        // Unlike `get_real_mass_vector`, retain unresolved or complex masses as unknown:
        // replacing either with zero could turn an inconclusive bound into a false rejection.
        let masses = graph.new_edgevec(|edge, _, _| {
            if matches!(edge.mass, EdgeMass::Zero) {
                Some(F(0.0))
            } else {
                edge.mass_value::<f64>(model, param_builder)
                    .and_then(|mass| IsZero::is_zero(&mass.im).then_some(mass.re))
            }
        });
        let Some(cut_boundary_mass) = Self::mass_sum(&self.cut_boundary_edges, &masses) else {
            return ThresholdCountertermStatus::PotentiallyExisting;
        };
        let Some(threshold_boundary_mass) = Self::mass_sum(&self.threshold_boundary_edges, &masses)
        else {
            return ThresholdCountertermStatus::PotentiallyExisting;
        };
        let e_cm = F(runtime_settings.kinematics.e_cm.abs());
        // The classified margins are differences of invariant masses squared, so the
        // dimensionless threshold is normalized by E_cm squared.
        let normalized_margin_tolerance = if normalized_margin_tolerance.is_finite() {
            normalized_margin_tolerance.abs()
        } else {
            DEFAULT_ESURFACE_EXISTENCE_THRESHOLD
        };
        let tolerance = F(normalized_margin_tolerance) * e_cm * e_cm;

        let maximum_cut_invariant = if self.cut_boundary_edges.len() == 1 {
            Some(cut_boundary_mass)
        } else {
            let incoming_invariant = Self::incoming_invariant_mass(graph, runtime_settings);
            let complement_mass = graph
                .iter_edges_of(&cut.cut)
                .map(|(_, edge_id, _)| edge_id)
                .filter(|edge_id| !self.cut_boundary_edges.contains(edge_id))
                .try_fold(F(0.0), |sum, edge_id| {
                    masses[edge_id].as_ref().map(|mass| sum + mass.abs())
                });
            incoming_invariant
                .zip(complement_mass)
                .and_then(|(invariant, complement_mass)| {
                    let maximum = invariant - complement_mass;
                    (maximum > F(0.0)).then_some(maximum)
                })
        };
        self.classify_from_invariant_bounds(
            &cut_boundary_mass,
            &threshold_boundary_mass,
            maximum_cut_invariant,
            &tolerance,
        )
    }
}

#[derive(Clone)]
struct TopologicalThresholdCandidate {
    topological_threshold_id: TopologicalThresholdId,
    left: SuBitGraph,
    cut: OrientedCut,
    right: SuBitGraph,
    threshold_esurface: Esurface,
    selected_orientation_esurface_id: Option<EsurfaceID>,
}

#[derive(Clone)]
struct ResolvedThresholdCountertermDraft {
    name: String,
    group_id: Option<usize>,
    origin: ThresholdCountertermOrigin,
    disable: bool,
    requested_subspace: Option<Vec<EdgeIndex>>,
    requested_parent_lmb: Option<Vec<EdgeIndex>>,
    multiplier: Option<ThresholdCountertermMultiplier>,
    subspace: SubspaceData,
    status: ThresholdCountertermStatus,
}

impl ResolvedThresholdCountertermDraft {
    fn is_compatible_with(
        &self,
        other: &Self,
        all_lmbs: &TiVec<LmbIndex, LoopMomentumBasis>,
    ) -> bool {
        self.name == other.name
            && self.group_id == other.group_id
            && self.disable == other.disable
            && self.multiplier == other.multiplier
            && self
                .subspace
                .has_equivalent_embedding(&other.subspace, all_lmbs)
    }
}

#[derive(Clone)]
struct PhysicalThresholdCountertermDraft {
    cut_id: CutId,
    cut_group_id: CutGroupId,
    cut_edges: Vec<EdgeIndex>,
    side: ThresholdCountertermSide,
    esurface_id: EsurfaceID,
    threshold_edges: Vec<EdgeIndex>,
    raised_threshold_id: RaisedEsurfaceId,
    raised_esurface_group: RaisedEsurfaceGroup,
    variants: Vec<ResolvedThresholdCountertermDraft>,
}

use derive_more::{From, Into};
#[derive(Clone, Encode, Decode)]
#[trait_decode(trait = GammaLoopContext)]
pub struct CrossSection {
    pub name: String,
    pub integrand: Option<ProcessIntegrand>,
    pub supergraphs: Vec<CrossSectionGraph>,
    pub external_particles: Vec<ArcParticle>,
    pub external_connections: Vec<ExternalConnection>,
    pub n_incmoming: usize,
    pub graph_group_structure: TiVec<GroupId, GraphGroup>,
}

impl CrossSection {
    pub fn plan_graph_group_selection(
        &self,
        spec: &GraphGroupSelectionSpec,
    ) -> Result<GraphGroupSelectionPlan> {
        if spec.mode() == GraphGroupSelectionMode::CrossSectionAmplitudeGraphs {
            return Err(eyre!(
                "`select --amplitude-graphs` for cross sections requires process and generation settings context."
            ));
        }
        if spec.has_raised_cut_rules() {
            return Err(eyre!(
                "Raised-cut signature selection for cross sections requires process and generation settings context."
            ));
        }
        spec.plan(&self.graph_group_structure, |graph_id| {
            self.supergraphs.get(graph_id).map(|graph| &graph.graph)
        })
    }

    pub fn plan_graph_group_selection_with_context(
        &self,
        spec: &GraphGroupSelectionSpec,
        model: &Model,
        process_definition: &ProcessDefinition,
        generation_settings: &GenerationSettings,
    ) -> Result<GraphGroupSelectionPlan> {
        match spec.mode() {
            GraphGroupSelectionMode::MasterGraphs => {
                spec.plan_with_analysis_contexts(
                    &self.graph_group_structure,
                    |graph_id| self.supergraphs.get(graph_id).map(|graph| &graph.graph),
                    |_master_graph_id, master_graph| {
                        Ok(vec![GraphSelectionSubject::whole_graph(master_graph)])
                    },
                    |master_graph_id, master_graph| {
                        self.cut_selection_subjects(
                            master_graph_id,
                            master_graph,
                            model,
                            process_definition,
                            generation_settings,
                        )
                    },
                    "Graph-group selection structural filters have no graph analysis subjects.",
                    "Raised-cut signature selection found no process-valid Cutkosky cuts to analyse.",
                )
            }
            GraphGroupSelectionMode::CrossSectionAmplitudeGraphs => {
                spec.plan_with_analysis_contexts(
                    &self.graph_group_structure,
                    |graph_id| self.supergraphs.get(graph_id).map(|graph| &graph.graph),
                    |master_graph_id, master_graph| {
                        let cross_section_graph =
                            self.supergraphs.get(master_graph_id).ok_or_else(|| {
                                eyre!(
                                    "Graph group refers to missing master supergraph id {}.",
                                    master_graph_id
                                )
                            })?;
                        let cuts = cross_section_graph.process_valid_cuts(
                            model,
                            process_definition,
                            generation_settings,
                        )?;
                        Ok(cuts
                            .into_iter()
                            .flat_map(|cut| cut.amplitude_side_subjects(master_graph))
                            .collect())
                    },
                    |master_graph_id, master_graph| {
                        self.cut_selection_subjects(
                            master_graph_id,
                            master_graph,
                            model,
                            process_definition,
                            generation_settings,
                        )
                    },
                    "`select --amplitude-graphs` found no process-valid Cutkosky cuts to analyse for structural filters.",
                    "Raised-cut signature selection found no process-valid Cutkosky cuts to analyse.",
                )
            }
        }
    }

    fn cut_selection_subjects<'a>(
        &'a self,
        master_graph_id: usize,
        master_graph: &'a Graph,
        model: &Model,
        process_definition: &ProcessDefinition,
        generation_settings: &GenerationSettings,
    ) -> Result<Vec<GraphCutSelectionSubject<'a>>> {
        let cross_section_graph = self.supergraphs.get(master_graph_id).ok_or_else(|| {
            eyre!(
                "Graph group refers to missing master supergraph id {}.",
                master_graph_id
            )
        })?;
        let cuts = cross_section_graph.process_valid_cuts(
            model,
            process_definition,
            generation_settings,
        )?;
        Ok(cuts
            .into_iter()
            .map(|cut| GraphCutSelectionSubject::new(master_graph, cut.cut.as_subgraph()))
            .collect())
    }

    pub fn amplitude_graph_signature_inventory(
        &self,
        model: &Model,
        process_definition: &ProcessDefinition,
        generation_settings: &GenerationSettings,
    ) -> Result<GraphSelectionSignatureInventory> {
        let subjects = self
            .graph_group_structure
            .iter()
            .flat_map(|group| group.into_iter().next())
            .map(|master_graph_id| {
                let cross_section_graph =
                    self.supergraphs.get(master_graph_id).ok_or_else(|| {
                        eyre!(
                            "Graph group refers to missing master supergraph id {}.",
                            master_graph_id
                        )
                    })?;
                let cuts = cross_section_graph.process_valid_cuts(
                    model,
                    process_definition,
                    generation_settings,
                )?;
                Ok(cuts
                    .into_iter()
                    .flat_map(|cut| cut.amplitude_side_subjects(&cross_section_graph.graph))
                    .collect::<Vec<_>>())
            })
            .collect::<Result<Vec<_>>>()?
            .into_iter()
            .flatten()
            .collect::<Vec<_>>();

        Ok(GraphSelectionSignatureInventory::from_analysis_subjects(
            subjects,
        ))
    }

    pub fn raised_cut_signature_inventory(
        &self,
        model: &Model,
        process_definition: &ProcessDefinition,
        generation_settings: &GenerationSettings,
    ) -> Result<RaisedCutSignatureInventory> {
        let subjects = self
            .graph_group_structure
            .iter()
            .flat_map(|group| group.into_iter().next())
            .map(|master_graph_id| {
                let cross_section_graph =
                    self.supergraphs.get(master_graph_id).ok_or_else(|| {
                        eyre!(
                            "Graph group refers to missing master supergraph id {}.",
                            master_graph_id
                        )
                    })?;
                let cuts = cross_section_graph.process_valid_cuts(
                    model,
                    process_definition,
                    generation_settings,
                )?;
                Ok(cuts
                    .into_iter()
                    .map(|cut| {
                        GraphCutSelectionSubject::new(
                            &cross_section_graph.graph,
                            cut.cut.as_subgraph(),
                        )
                    })
                    .collect::<Vec<_>>())
            })
            .collect::<Result<Vec<_>>>()?
            .into_iter()
            .flatten()
            .collect::<Vec<_>>();

        Ok(RaisedCutSignatureInventory::from_cut_subjects(subjects))
    }

    pub fn validate_graph_group_selection_plan(
        &self,
        plan: &GraphGroupSelectionPlan,
    ) -> Result<()> {
        for &old_group_id in plan.retained_group_ids() {
            plan.new_group_id_for_old(old_group_id).ok_or_else(|| {
                eyre!(
                    "Selection plan is missing compact group id for old group {}.",
                    old_group_id.0
                )
            })?;
            if old_group_id.0 >= self.graph_group_structure.len() {
                return Err(eyre!(
                    "Selection plan refers to missing graph group {}.",
                    old_group_id.0
                ));
            }
            let group = &self.graph_group_structure[old_group_id];
            for old_graph_id in group {
                if old_graph_id >= self.supergraphs.len() {
                    return Err(eyre!(
                        "Graph group {} refers to missing graph id {}.",
                        old_group_id.0,
                        old_graph_id
                    ));
                }
            }
            let master = group.master();
            if master >= self.supergraphs.len() {
                return Err(eyre!(
                    "Graph group {} refers to missing master graph id {}.",
                    old_group_id.0,
                    master
                ));
            }
        }

        Ok(())
    }

    pub fn apply_graph_group_selection(&mut self, plan: &GraphGroupSelectionPlan) -> Result<()> {
        self.validate_graph_group_selection_plan(plan)?;

        let mut old_graph_to_new_group = vec![None; self.supergraphs.len()];
        let mut old_graph_is_master = vec![false; self.supergraphs.len()];
        for &old_group_id in plan.retained_group_ids() {
            let new_group_id = plan.new_group_id_for_old(old_group_id).ok_or_else(|| {
                eyre!(
                    "Selection plan is missing compact group id for old group {}.",
                    old_group_id.0
                )
            })?;
            let group = &self.graph_group_structure[old_group_id];
            for old_graph_id in group {
                old_graph_to_new_group[old_graph_id] = Some(new_group_id);
            }
            old_graph_is_master[group.master()] = true;
        }

        let mut new_supergraphs = self
            .supergraphs
            .iter()
            .cloned()
            .enumerate()
            .filter_map(|(old_graph_id, mut graph)| {
                let new_group_id = old_graph_to_new_group[old_graph_id]?;
                graph.graph.group_id = Some(new_group_id);
                graph.graph.is_group_master = old_graph_is_master[old_graph_id];
                if let Some(multi_channeling_setup) = &mut graph.derived_data.multi_channeling_setup
                {
                    multi_channeling_setup.graph.group_id = Some(new_group_id);
                    multi_channeling_setup.graph.is_group_master =
                        old_graph_is_master[old_graph_id];
                }
                Some(graph)
            })
            .collect::<Vec<_>>();

        let mut parsed_graphs = new_supergraphs
            .iter()
            .map(|graph| graph.graph.clone())
            .collect::<Vec<_>>();
        let new_graph_group_structure = complete_group_parsing(&mut parsed_graphs)?;
        for (graph, parsed_graph) in new_supergraphs.iter_mut().zip(parsed_graphs) {
            graph.graph.group_id = parsed_graph.group_id;
            graph.graph.is_group_master = parsed_graph.is_group_master;
            if let Some(multi_channeling_setup) = &mut graph.derived_data.multi_channeling_setup {
                multi_channeling_setup.graph.group_id = graph.graph.group_id;
                multi_channeling_setup.graph.is_group_master = graph.graph.is_group_master;
            }
        }

        self.supergraphs = new_supergraphs;
        self.graph_group_structure = new_graph_group_structure;
        self.integrand = None;
        Ok(())
    }

    fn storage_path(&self, base: &Path) -> PathBuf {
        base.join(&self.name)
    }

    pub fn export_standalone(
        &self,
        path: impl AsRef<Path>,
        settings: &crate::processes::StandaloneExportSettings,
    ) -> Result<()> {
        if let Some(integrand) = &self.integrand {
            integrand.export_standalone(path, settings)?;
        } else {
            return Err(eyre!(
                "Cannot export standalone cross section {} without integrand",
                self.name
            ));
        }

        Ok(())
    }

    #[allow(dead_code)]
    pub(crate) fn write_dot<W: std::io::Write>(
        &self,
        writer: &mut W,
        settings: &DotExportSettings,
    ) -> Result<(), std::io::Error> {
        for graph in &self.supergraphs {
            graph.write_dot(writer, settings)?;
            writeln!(writer)?;
        }
        Ok(())
    }

    pub fn write_dot_fmt<W: std::fmt::Write>(
        &self,
        writer: &mut W,
        settings: &DotExportSettings,
    ) -> Result<(), std::fmt::Error> {
        for graph in &self.supergraphs {
            graph.write_dot_fmt(writer, settings)?;
            writeln!(writer)?;
        }
        Ok(())
    }

    pub(crate) fn new(name: String) -> Self {
        Self {
            name,
            integrand: None,
            supergraphs: vec![],
            external_connections: vec![],
            external_particles: vec![],
            n_incmoming: 0,
            graph_group_structure: TiVec::new(),
        }
    }

    pub fn from_graph_list(name: String, mut graphs: Vec<Graph>, _model: &Model) -> Result<Self> {
        let mut cross_section = CrossSection::new(name);
        cross_section.graph_group_structure = complete_group_parsing(&mut graphs)?;
        // println!("group structure: {:?}", cross_section.graph_group_structure);

        for group in &cross_section.graph_group_structure {
            for graph_id in group.into_iter().skip(1) {
                let graph = &graphs[graph_id];
                if !graph.threshold_counterterms.autogenerated {
                    return Err(eyre!(
                        "Cross-section graph '{}' is a non-master member of graph group '{}', but declares threshold_counterterms metadata. Threshold metadata on non-master cross-section graphs is unsupported: cross-section threshold solve groups are graph-local. Put the metadata on the master graph or import this graph in its own group.",
                        graph.name,
                        graphs[group.master()].name,
                    ));
                }
            }
        }

        for cross_section_graph in graphs {
            // cross_section_graph.param_builder =  ParamBuilder::new(&cross_section_graph, model);
            cross_section.add_supergraph(cross_section_graph)?;
        }
        Ok(cross_section)
    }

    pub(crate) fn warm_up(&mut self, model: &Model) -> Result<()> {
        if let Some(integrand) = &mut self.integrand {
            integrand.warm_up(model)
        } else {
            Err(eyre!(
                "Cannot warm up amplitude {} without integrand",
                self.name
            ))
        }
    }

    fn add_supergraph(&mut self, supergraph: Graph) -> Result<()> {
        if self.external_particles.is_empty() {
            let external_particles = supergraph.get_external_partcles();
            if !external_particles.len().is_multiple_of(2) {
                return Err(eyre!(
                    "expected even number of externals for forward scattering graph"
                ));
            }
            self.external_particles = external_particles;
            self.n_incmoming = self.external_particles.len() / 2;
        } else if self.external_particles != supergraph.get_external_partcles() {
            return Err(eyre!(
                "attempt to add supergraph with differnt external particles"
            ));
        }

        let cross_section_graph = CrossSectionGraph::new(supergraph);
        self.supergraphs.push(cross_section_graph);

        // TODO: validate that the graph is compatible
        Ok(())
    }

    pub fn preprocess(
        &mut self,
        model: &Model,
        process_definition: &ProcessDefinition,
        global_settings: &GenerationSettings,
        runtime_default: LockedRuntimeSettings,
        generation_pool: &ThreadPool,
    ) -> Result<Vec<NamedGraphGenerationReport>> {
        global_settings.validate_for_process(GenerationType::CrossSection)?;
        let integrand_name = self.name.clone();
        generation_progress::begin_phase(
            GenerationProgressPhase::GraphPreprocessing,
            GenerationProcessKind::CrossSection,
            &process_definition.folder_name,
            &integrand_name,
            self.supergraphs.len(),
            None,
        );
        generation_pool.install(|| {
            self.supergraphs
                .par_iter_mut()
                .map(|supergraph| {
                    if crate::is_interrupted() {
                        return Err(eyre!("Generation interrupted by user"));
                    }
                    let graph_name = supergraph.graph.name.clone();
                    generation_progress::graph_started(
                        GenerationProcessKind::CrossSection,
                        &integrand_name,
                        &graph_name,
                        None,
                    );
                    let stats = supergraph.preprocess(
                        model,
                        process_definition,
                        global_settings,
                        runtime_default,
                    )?;
                    if crate::is_interrupted() {
                        return Err(eyre!("Generation interrupted by user"));
                    }
                    generation_progress::graph_finished(
                        GenerationProcessKind::CrossSection,
                        &integrand_name,
                        &graph_name,
                        &stats,
                        None,
                    );
                    Ok(NamedGraphGenerationReport {
                        integrand_name: integrand_name.clone(),
                        graph_name,
                        stats,
                    })
                })
                .collect::<Result<Vec<_>>>()
        })
    }

    pub fn build_integrand(
        &mut self,
        model: &Model,
        process_definition: &ProcessDefinition,
        global_settings: &GlobalSettings,
        runtime_default: LockedRuntimeSettings,
        generation_pool: &ThreadPool,
    ) -> Result<Vec<NamedGraphGenerationReport>> {
        let process_name = process_definition.folder_name.as_str();
        let started = std::time::Instant::now();
        crate::debug_tags!(#generation, #profile, #graph, #summary;
            stage = "cross_section_build_integrand_start",
            integrand = %self.name,
            graph_count = self.supergraphs.len(),
            "Generation timing milestone"
        );
        if crate::is_interrupted() {
            return Err(eyre!("Generation interrupted by user"));
        }
        let integrand_name = self.name.clone();
        let total_cuts = self
            .supergraphs
            .iter()
            .map(|sg| {
                sg.derived_data
                    .cut_group_data
                    .cut_groups
                    .iter()
                    .map(|cut_group| cut_group.cuts.len())
                    .sum::<usize>()
            })
            .sum();
        generation_progress::begin_phase(
            GenerationProgressPhase::GraphGeneration,
            GenerationProcessKind::CrossSection,
            process_name,
            &integrand_name,
            self.supergraphs.len(),
            Some(total_cuts),
        );
        let mut graph_reports = Vec::new();
        let terms = generation_pool.install(|| {
            self.supergraphs
                .par_iter_mut()
                .map(|sg| {
                    if crate::is_interrupted() {
                        return Err(eyre!("Generation interrupted by user"));
                    }
                    let graph_started = std::time::Instant::now();
                    crate::debug_tags!(#generation, #profile, #graph, #summary;
                        stage = "generate_term_for_graph_start",
                        integrand = %integrand_name,
                        graph = %sg.graph.name,
                        "Generation timing milestone"
                    );
                    generation_progress::graph_started(
                        GenerationProcessKind::CrossSection,
                        &integrand_name,
                        &sg.graph.name,
                        Some(sg.cuts.len()),
                    );
                    let _progress_context_guard = generation_progress::enter_progress_context(
                        format!("{} / {} cuts", sg.graph.name, sg.cuts.len()),
                    );
                    let (term, mut stats) = sg.generate_term_for_graph(model, global_settings)?;
                    if crate::is_interrupted() {
                        return Err(eyre!("Generation interrupted by user"));
                    }
                    stats.total_time += graph_started.elapsed();
                    crate::debug_tags!(#generation, #profile, #graph, #summary;
                        stage = "generate_term_for_graph_done",
                        integrand = %integrand_name,
                        graph = %sg.graph.name,
                        elapsed_ms = graph_started.elapsed().as_secs_f64() * 1000.0,
                        "Generation timing milestone"
                    );
                    generation_progress::graph_finished(
                        GenerationProcessKind::CrossSection,
                        &integrand_name,
                        &sg.graph.name,
                        &stats,
                        None,
                    );
                    Ok((
                        term,
                        NamedGraphGenerationReport {
                            integrand_name: integrand_name.clone(),
                            graph_name: sg.graph.name.clone(),
                            stats,
                        },
                    ))
                })
                .collect::<Result<Vec<_>>>()
        })?;
        if crate::is_interrupted() {
            return Err(eyre!("Generation interrupted by user"));
        }
        for (_, report) in &terms {
            graph_reports.push(report.clone());
        }
        let mut terms = terms.into_iter().map(|(term, _)| term).collect::<Vec<_>>();

        for group in self.graph_group_structure.iter() {
            let master = group.master();
            let mc_of_master = self.supergraphs[master]
                .derived_data
                .multi_channeling_setup
                .as_ref()
                .unwrap();

            for graph_id in group.into_iter() {
                terms[graph_id].multi_channeling_setup = mc_of_master.clone();
            }
        }

        let backend_started = std::time::Instant::now();
        let graph_count = terms.len();
        crate::debug_tags!(#generation, #profile, #compile, #graph, #summary;
            stage = "prepare_runtime_backends_start",
            integrand = %self.name,
            graph_count,
            elapsed_ms = started.elapsed().as_secs_f64() * 1000.0,
            "Generation timing milestone"
        );
        generation_progress::backend_started(
            GenerationProcessKind::CrossSection,
            &self.name,
            graph_count,
        );
        let mut cross_section_integrand = CrossSectionIntegrand {
            settings: runtime_default.into(),
            data: CrossSectionIntegrandData {
                compilation: global_settings
                    .generation
                    .compile
                    .frozen_mode(&global_settings.generation.evaluator),
                loop_cache_id: 0,
                external_cache_id: 0,
                base_external_cache_id: 0,
                rotations: None,
                symmetrize_left_right_states: process_definition.symmetrize_left_right_states,
                name: self.name.clone(),
                external_connections: self.external_connections.clone(),
                n_incoming: self.n_incmoming,
                // polarizations,
                graph_terms: terms,
                graph_group_structure: self.graph_group_structure.clone(),
                graph_to_group_id: graph_to_group_id_for_group_structure(
                    &self.graph_group_structure,
                ),
                explicit_orientation_sum_only: global_settings
                    .generation
                    .explicit_orientation_sum_only,
            },
            event_processing_runtime: Default::default(),
            active_f64_backend: Default::default(),
        };
        let compile_times = cross_section_integrand
            .prepare_runtime_backends_after_generation_with_compile_times()?;
        crate::debug_tags!(#generation, #profile, #compile, #graph, #summary;
            stage = "prepare_runtime_backends_done",
            integrand = %self.name,
            graph_count,
            elapsed_ms = backend_started.elapsed().as_secs_f64() * 1000.0,
            total_elapsed_ms = started.elapsed().as_secs_f64() * 1000.0,
            "Generation timing milestone"
        );
        generation_progress::backend_finished(
            GenerationProcessKind::CrossSection,
            &self.name,
            backend_started.elapsed(),
        );
        for (report, compile_time) in graph_reports.iter_mut().zip(compile_times) {
            report.stats.evaluator_compile_time += compile_time;
            report.stats.total_time += compile_time;
        }

        self.integrand = Some(ProcessIntegrand::CrossSection(cross_section_integrand));
        crate::debug_tags!(#generation, #profile, #graph, #summary;
            stage = "cross_section_build_integrand_done",
            integrand = %self.name,
            graph_count = self.supergraphs.len(),
            elapsed_ms = started.elapsed().as_secs_f64() * 1000.0,
            "Generation timing milestone"
        );
        Ok(graph_reports)
    }

    pub fn compile(
        &mut self,
        path: impl AsRef<Path>,
        override_existing: bool,
        thread_pool: &ThreadPool,
    ) -> Result<Vec<NamedGraphGenerationReport>> {
        info!("Compiling cross section {}", self.name);
        let p = self.storage_path(path.as_ref());

        let result = fs::create_dir_all(&p).with_context(|| {
            format!(
                "Trying to create directory to compile cross section {}",
                p.display()
            )
        });

        if override_existing {
            result?;
        }

        if let Some(integrand) = &mut self.integrand {
            let compile_times = integrand.compile(&p, override_existing, thread_pool)?;
            return Ok(compile_times
                .into_iter()
                .map(|(graph_name, duration)| NamedGraphGenerationReport {
                    integrand_name: self.name.clone(),
                    graph_name,
                    stats: GraphGenerationStats {
                        total_time: duration,
                        evaluator_compile_time: duration,
                        ..GraphGenerationStats::default()
                    },
                })
                .collect());
        }
        Ok(Vec::new())
    }

    pub fn save(&mut self, path: impl AsRef<Path>, override_existing: bool) -> Result<()> {
        let p = self.storage_path(path.as_ref());
        let r = fs::create_dir_all(&p).with_context(|| {
            format!(
                "Trying to create directory to save cross section {}",
                p.display()
            )
        });

        if override_existing {
            r?;
        }

        let integrand = self.integrand.take();
        if let Some(integrand) = &integrand {
            integrand.save(&p, override_existing)?;
        }

        let binary = bincode::encode_to_vec(&(*self), bincode::config::standard())?;
        if override_existing {
            fs::write(p.join("cs.bin"), &binary)?;
        } else {
            let mut file = File::create_new(p.join("cs.bin"))?;
            file.write_all(&binary)?;
        }

        self.integrand = integrand;
        Ok(())
    }

    pub(crate) fn load(path: impl AsRef<Path>, context: GammaLoopContextContainer) -> Result<Self> {
        let binary = fs::read(path.as_ref().join("cs.bin"))?;
        let (mut cs, _): (Self, _) =
            bincode::decode_from_slice_with_context(&binary, bincode::config::standard(), context)?;

        if path.as_ref().join("integrand").exists() {
            let integrand = CrossSectionIntegrand::load(path.as_ref().join("integrand"), context)?;
            cs.integrand = Some(ProcessIntegrand::CrossSection(integrand));
        }

        Ok(cs)
    }
}

#[derive(Clone, bincode::Encode, bincode::Decode)]
pub struct CrossSectionCut {
    pub cut: OrientedCut,
    pub left: SuBitGraph,
    pub right: SuBitGraph,
}

impl CrossSectionCut {
    fn amplitude_side_subjects<'a>(&self, graph: &'a Graph) -> [GraphSelectionSubject<'a>; 2] {
        let cut_edges = self.cut.as_subgraph();
        [
            GraphSelectionSubject::cut_side_amplitude_subgraph(
                graph,
                self.left.subtract(&cut_edges),
            ),
            GraphSelectionSubject::cut_side_amplitude_subgraph(
                graph,
                self.right.subtract(&cut_edges),
            ),
        ]
    }

    pub(crate) fn is_s_channel(&self, cross_section_graph: &CrossSectionGraph) -> Result<bool> {
        let nodes_of_left_cut: Vec<_> = cross_section_graph
            .graph
            .underlying
            .iter_nodes_of(&self.left)
            .map(|(nid, _, _)| nid)
            .collect();

        let left_node = cross_section_graph
            .graph
            .underlying
            .combine_to_single_hedgenode(&nodes_of_left_cut);
        let res = left_node.includes(&cross_section_graph.source_nodes)
            && cross_section_graph.target_nodes.weakly_disjoint(&left_node);

        if !res {
            warn!("s channel check wrong");
        }

        Ok(true)
    }

    pub(crate) fn is_valid_for_process(
        &self,
        cross_section_graph: &CrossSectionGraph,
        process: &ProcessDefinition,
        model: &Model,
    ) -> Result<bool> {
        if self.is_s_channel(cross_section_graph)? {
            let cut_content_builder = self
                .cut
                .iter_edges(&cross_section_graph.graph.underlying)
                .filter_map(|(orientation, edge_data)| {
                    Some(if orientation == Orientation::Reversed {
                        edge_data.data.particle()?.get_anti_particle(model)
                    } else {
                        edge_data.data.particle()?.clone()
                    })
                })
                .collect_vec();

            let covariant_states = process.covariant_cut_states(model)?;
            let any_pdg_list_passes = covariant_states
                .iter()
                .map(|x| {
                    x.iter()
                        .map(|pdg| model.get_particle_from_pdg(*pdg as isize))
                })
                .any(|particle_content| {
                    let mut cut_content = cut_content_builder.clone();
                    debug!(
                        "cut content: {:?}",
                        cut_content.iter().map(|p| p.name.clone()).collect_vec()
                    );

                    for particle in particle_content {
                        if let Some(index) = cut_content.iter().position(|p| p == &particle) {
                            cut_content.remove(index);
                        } else {
                            debug!("wrong particles");
                            return false;
                        }
                    }

                    let (n_unresolved, unresolved_cut_content) =
                        process.unresolved_cut_content(model);

                    if cut_content.len() > n_unresolved {
                        debug!(" too many unresolved particles");
                        return false;
                    }

                    if !cut_content
                        .iter()
                        .all(|particle| unresolved_cut_content.contains(particle))
                    {
                        debug!("wrong unresolved particles");
                        return false;
                    }

                    true
                });

            if !any_pdg_list_passes {
                debug!("wrong pdg list");
                return Ok(false);
            }

            let amplitude_couplings = process.amplitude_filters.get_coupling_orders();
            let amplitude_loop_count = process.amplitude_filters.get_loop_count_range();

            if let Some((min_loop, max_loop)) = amplitude_loop_count {
                let loop_range = min_loop..=max_loop;
                let left_internal_subgraph = InternalSubGraph::cleaned_filter_pessimist(
                    self.left.clone(),
                    &cross_section_graph.graph.underlying,
                );

                let right_internal_subgraph = InternalSubGraph::cleaned_filter_pessimist(
                    self.right.clone(),
                    &cross_section_graph.graph.underlying,
                );

                let left_loop = cross_section_graph
                    .graph
                    .underlying
                    .cyclotomatic_number(&left_internal_subgraph);

                let right_loop = cross_section_graph
                    .graph
                    .underlying
                    .cyclotomatic_number(&right_internal_subgraph);

                let total_loops = left_loop + right_loop;

                if !loop_range.contains(&total_loops) {
                    debug!("incorrect loop count");
                    return Ok(false);
                }
            }

            if amplitude_couplings.is_some() {
                todo!("waiting for update")
            }

            Ok(true)
        } else {
            debug!("cut is not s channel");
            Ok(false)
        }
    }
}

#[derive(
    Debug, Clone, Serialize, Decode, Deserialize, From, Into, Hash, PartialEq, Copy, Eq, Encode,
)]
pub struct CutId(pub usize);

impl Display for CutId {
    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
        write!(f, "{}", self.0)
    }
}

#[derive(Clone, Encode, Decode)]
#[trait_decode(trait = GammaLoopContext)]
pub struct CrossSectionGraph {
    pub graph: Graph,
    pub source_nodes: HedgeNode,
    pub target_nodes: HedgeNode,
    pub cuts: TiVec<CutId, CrossSectionCut>,
    pub cut_esurface: TiVec<CutId, Esurface>,
    pub cut_esurface_id_map: TiVec<CutId, EsurfaceID>,
    pub derived_data: CrossSectionDerivedData,
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct CutkoskyCutCount {
    pub candidate_st_cuts: usize,
    pub multi_edge_candidate_cuts: usize,
    pub process_compatible_cuts: usize,
    pub selected_cuts: usize,
}

impl CrossSectionGraph {
    pub(crate) fn new(graph: Graph) -> Self {
        let (source_node, target_node) = graph.get_source_and_target();

        Self {
            graph,
            source_nodes: source_node,
            target_nodes: target_node,
            cuts: TiVec::new(),
            cut_esurface: TiVec::new(),
            cut_esurface_id_map: TiVec::new(),
            derived_data: CrossSectionDerivedData::new_empty(),
        }
    }

    pub(crate) fn apply_spin_sum(
        &mut self,
        model: &Model,
        generation_settings: &GenerationSettings,
        locked_runtime_settings: &LockedRuntimeSettings,
    ) -> Result<()> {
        for (extid, hel) in locked_runtime_settings.helicities().iter().enumerate() {
            //println!("{extid},{hel:?}");
            let eid = self.graph.loop_momentum_basis.ext_edges[ExternalIndex(extid)];

            let Some(p) = self.graph.underlying[eid].particle() else {
                continue;
            };

            match hel {
                Helicity::Summed => {
                    let Some(p) = p.polarization_sum(
                        eid,
                        false,
                        generation_settings.vector_polarization_sum_gauge,
                    )?
                    else {
                        continue;
                    };
                    self.graph.global_prefactor.projector =
                        self.graph.global_prefactor.projector.replace_multiple(&[p]);
                }
                Helicity::SummedAveraged => {
                    let Some(p) = p.polarization_sum(
                        eid,
                        true,
                        generation_settings.vector_polarization_sum_gauge,
                    )?
                    else {
                        continue;
                    };
                    self.graph.global_prefactor.projector =
                        self.graph.global_prefactor.projector.replace_multiple(&[p]);
                }
                _ => {}
            }
        }

        self.graph.polarizations = self.graph.global_prefactor.polarizations();
        self.graph.param_builder = ParamBuilder::new(
            &self.graph,
            model,
            &self.graph.loop_momentum_basis,
            &self.graph.param_builder.pairs.additional_params.params,
        )?;
        Ok(())
    }

    pub(crate) fn preprocess(
        &mut self,
        model: &Model,
        process_definition: &ProcessDefinition,
        settings: &GenerationSettings,
        runtime_default: LockedRuntimeSettings,
    ) -> Result<GraphGenerationStats> {
        let preprocess_started = std::time::Instant::now();
        let mut stats = GraphGenerationStats::default();
        self.graph.validate_real_masses(model)?;
        process_definition.covariant_cut_states(model)?;
        self.derived_data.covariant_cut_representatives =
            process_definition.covariant_cut_representatives(model);
        self.apply_spin_sum(model, settings, &runtime_default)?;
        debug_tags!(#generation; "generating cuts");
        self.generate_cuts(model, process_definition, settings)?;
        debug_tags!(#generation; "generating esurfaces corresponding to cuts");
        self.generate_esurface_cuts();
        debug_tags!(#generation; "generating cff");
        stats.merge_in_place(&self.generate_cff(settings)?);
        debug_tags!(#generation; "building lmbs");
        self.build_lmbs()?;
        debug_tags!(#generation; "building multi channeling channels");

        if self.graph.is_group_master {
            self.build_multi_channeling_channels(settings.override_lmb_heuristics)?;
        }

        let vk = crate::utils::vakint()?;
        debug_tags!(#generation; "building parametric integrand");
        let cff_energy_degree_bound_reports = Mutex::new(Vec::new());
        self.build_parametric_integrand(settings, vk, &cff_energy_degree_bound_reports)?;
        //self.build_parametric_integrand_cut_groups(settings)?;

        let threshold_candidates = self.topological_threshold_candidates(settings)?;
        self.derived_data.threshold_candidate_esurface_ids = threshold_candidates
            .iter()
            .filter_map(|candidate| candidate.selected_orientation_esurface_id)
            .sorted()
            .dedup()
            .collect();
        self.derived_data.resolved_threshold_counterterms = None;
        self.derived_data.threshold_counterterms = TiVec::new();
        self.derived_data.subspace_data = TiVec::new();

        self.derived_data.cut_threshold_associations =
            ti_vec![CutThresholdCountertermAssociations::default(); self.cuts.len()];

        if settings.threshold_subtraction.enable_thresholds {
            debug_tags!(#generation, #subtraction; "building threshold counterterm");
            self.build_subspace_data()?;
            let runtime_settings: RuntimeSettings = runtime_default.into();
            if let Err(error) = self.build_threshold_counterterm(
                model,
                settings,
                &runtime_settings,
                &threshold_candidates,
                vk,
                &cff_energy_degree_bound_reports,
            ) {
                self.derived_data.resolved_threshold_counterterms = None;
                self.derived_data.threshold_counterterms = TiVec::new();
                self.derived_data.subspace_data = TiVec::new();
                self.derived_data.cut_threshold_associations =
                    ti_vec![CutThresholdCountertermAssociations::default(); self.cuts.len()];
                return Err(error);
            }
        }

        stats.cff_energy_degree_bound_reports = cff_energy_degree_bound_reports
            .into_inner()
            .unwrap_or_else(|poisoned| poisoned.into_inner());

        stats.total_time += preprocess_started.elapsed();
        Ok(stats)
    }

    #[allow(dead_code)]
    pub(crate) fn write_dot<W: std::io::Write>(
        &self,
        writer: &mut W,
        settings: &DotExportSettings,
    ) -> Result<(), std::io::Error> {
        if let Some(graph) = self.graph_with_materialized_threshold_counterterms(settings) {
            graph.dot_serialize_io(writer, settings)
        } else {
            self.graph.dot_serialize_io(writer, settings)
        }
    }

    pub(crate) fn write_dot_fmt<W: std::fmt::Write>(
        &self,
        writer: &mut W,
        settings: &DotExportSettings,
    ) -> Result<(), std::fmt::Error> {
        if let Some(graph) = self.graph_with_materialized_threshold_counterterms(settings) {
            graph.dot_serialize_fmt(writer, settings)
        } else {
            self.graph.dot_serialize_fmt(writer, settings)
        }
    }

    fn graph_with_materialized_threshold_counterterms(
        &self,
        settings: &DotExportSettings,
    ) -> Option<Graph> {
        // Non-master cross-section members cannot accept threshold metadata on reimport.
        if !settings.include_autogenerated_fields
            || (self.graph.group_id.is_some() && !self.graph.is_group_master)
        {
            return None;
        }
        let resolved = self.derived_data.resolved_threshold_counterterms.as_ref()?;
        let all_lmbs = self.derived_data.lmbs.as_ref()?;
        let mut graph = self.graph.clone();
        graph.threshold_counterterms = crate::graph::autogen::Autogen::explicit(
            resolved.materialized_spec(&self.graph.threshold_counterterms, all_lmbs),
        );
        Some(graph)
    }

    fn generate_cff(&mut self, settings: &GenerationSettings) -> Result<GraphGenerationStats> {
        self.graph.set_medium_mode(settings.medium.mode)?;
        self.graph.ensure_energy_convergent_cycles(
            &self
                .graph
                .no_dummy()
                .subtract(&self.graph.initial_state_cut),
        )?;
        let canonize_esurface = self
            .graph
            .get_esurface_canonization(&self.graph.loop_momentum_basis);

        let contract_edges = self
            .graph
            .iter_edges_of(
                &self
                    .graph
                    .tree_edges
                    .subtract(&self.graph.initial_state_cut),
            )
            .map(|x| x.1)
            .collect_vec();

        // The options provision maps from the complete factorized numerator.
        // The numerator itself remains caller-owned and is mapped only when
        // the corresponding orientation integrand is assembled.
        let options = self.graph.production_cff_3d_expression_options(settings)?;
        let mut global_cff = self.graph.generate_3d_expression_for_integrand(
            &contract_edges,
            &canonize_esurface,
            &options,
            None,
        )?;
        let raised_edge_groups = self.graph.get_raised_edge_groups();
        normalize_three_d_expression_cut_support_with_raised_edge_groups(
            &mut global_cff.expression,
            &raised_edge_groups,
        );

        let generated_esurfaces = &global_cff.expression.surfaces.esurface_cache;
        let normalized_generated_esurfaces = generated_esurfaces
            .iter()
            .map(|esurface| {
                Graph::normalize_esurface_with_raised_edge_groups(esurface, &raised_edge_groups)
            })
            .collect::<TiVec<EsurfaceID, _>>();
        let cut_esurface_map = self
            .cut_esurface
            .iter()
            .map(|esurface| {
                let normalized =
                    Graph::normalize_esurface_with_raised_edge_groups(esurface, &raised_edge_groups);
                generated_esurfaces
                    .iter()
                    .position(|candidate| candidate == esurface)
                    .or_else(|| {
                        normalized_generated_esurfaces
                            .iter()
                            .position(|candidate| candidate == &normalized)
                    })
                    .map(EsurfaceID)
                    .ok_or_else(|| {
                        eyre!(
                            "Cutkosky-cut E-surface {esurface:?} for graph {} is absent from its generated CFF expression catalogue (normalized: {normalized:?})",
                            self.graph.name
                        )
                    })
            })
            .collect::<Result<_>>()?;

        self.cut_esurface_id_map = cut_esurface_map;

        let esurface_raised_data = self
            .graph
            .determine_raised_esurfaces_from_expression(&global_cff.expression);

        let (cut_group_data, cut_group_stats) = CutGroupData::new_from_esurface(
            &esurface_raised_data,
            &self.cut_esurface_id_map,
            &settings.evaluator,
        );

        self.derived_data.global_cff_expression = Some(global_cff);
        self.derived_data.cut_group_data = cut_group_data;

        Ok(cut_group_stats)
    }

    pub(crate) fn process_valid_cuts(
        &self,
        model: &Model,
        process_definition: &ProcessDefinition,
        settings: &GenerationSettings,
    ) -> Result<TiVec<CutId, CrossSectionCut>> {
        if !self.cuts.is_empty() {
            return Ok(self.cuts.clone());
        }
        self.compute_process_valid_cuts(model, process_definition, settings)
            .map(|(_, cuts)| cuts)
    }

    pub fn cutkosky_cut_count_for_process(
        &self,
        model: &Model,
        process_definition: &ProcessDefinition,
        settings: &GenerationSettings,
    ) -> Result<CutkoskyCutCount> {
        self.compute_process_valid_cuts(model, process_definition, settings)
            .map(|(cut_count, _)| cut_count)
    }

    fn compute_process_valid_cuts(
        &self,
        model: &Model,
        process_definition: &ProcessDefinition,
        settings: &GenerationSettings,
    ) -> Result<(CutkoskyCutCount, TiVec<CutId, CrossSectionCut>)> {
        let all_st_cuts = self.graph.all_st_cuts_for_cs(
            self.source_nodes.clone(),
            self.target_nodes.clone(),
            &self.graph.get_initial_state_tree(),
        );
        let candidate_st_cuts = all_st_cuts.len();

        let multi_edge_cuts: TiVec<CutId, CrossSectionCut> = all_st_cuts
            .into_iter()
            .map(|(left, cut, right)| CrossSectionCut { cut, left, right })
            .filter(|cut| cut.cut.nedges(&self.graph) > 1)
            .collect();
        let multi_edge_candidate_cuts = multi_edge_cuts.len();

        let mut cuts: TiVec<CutId, CrossSectionCut> = multi_edge_cuts
            .into_iter()
            .filter_map(
                |cut| match cut.is_valid_for_process(self, process_definition, model) {
                    Ok(true) => Some(Ok(cut)),
                    Ok(false) => None,
                    Err(e) => Some(Err(e)),
                },
            )
            .collect::<Result<_>>()?;
        let process_compatible_cuts = cuts.len();

        cuts.sort_by(|a, b| a.cut.cmp(&b.cut));

        if !settings.force_cuts.is_empty() {
            let force_cuts_sorted = settings
                .force_cuts
                .clone()
                .into_iter()
                .map(|cut_edges| cut_edges.into_iter().sorted().collect_vec())
                .collect_vec();

            cuts.retain(|cut| {
                let edges_in_cut = self
                    .graph
                    .iter_edges_of(&cut.cut)
                    .map(|(_, _, e)| e.data.name.value.clone())
                    .sorted()
                    .collect_vec();

                force_cuts_sorted.contains(&edges_in_cut)
            });
        }

        let cut_count = CutkoskyCutCount {
            candidate_st_cuts,
            multi_edge_candidate_cuts,
            process_compatible_cuts,
            selected_cuts: cuts.len(),
        };
        Ok((cut_count, cuts))
    }

    fn generate_cuts(
        &mut self,
        model: &Model,
        process_definition: &ProcessDefinition,
        settings: &GenerationSettings,
    ) -> Result<()> {
        debug_tags!(#generation, #profile, #graph;
            stage = "cross_section_generate_cuts_start",
            graph = %self.graph.name,
            "Cut discovery timing milestone"
        );
        let started = std::time::Instant::now();
        let (cut_count, cuts) =
            self.compute_process_valid_cuts(model, process_definition, settings)?;
        self.cuts = cuts;
        debug_tags!(#generation, #profile, #graph;
            stage = "cross_section_generate_cuts_done",
            graph = %self.graph.name,
            st_cut_count = cut_count.candidate_st_cuts,
            valid_cut_count = self.cuts.len(),
            elapsed_ms = started.elapsed().as_secs_f64() * 1000.0,
            "Cut discovery timing milestone"
        );
        generation_progress::cuts_discovered(
            "",
            &self.graph.name,
            cut_count.candidate_st_cuts,
            self.cuts.len(),
        );

        Ok(())
    }

    fn generate_esurface_cuts(&mut self) {
        debug!("generating esurfaces for cuts");

        let esurfaces: TiVec<CutId, Esurface> = self
            .cuts
            .iter()
            .map(|cut| {
                Esurface::new_from_cut_left(
                    &self.graph.underlying,
                    cut,
                    Some(&self.graph.initial_state_cut),
                )
            })
            .collect();

        debug!("generated esurfaces {:?}", esurfaces);

        self.cut_esurface = esurfaces;
    }

    pub(crate) fn build_parametric_integrand(
        &mut self,
        settings: &GenerationSettings,
        vakint: &Vakint,
        cff_energy_degree_bound_reports: &Mutex<Vec<CffEnergyDegreeBoundReport>>,
    ) -> Result<()> {
        self.derived_data.cut_paramatric_integrand =
            self.build_integrand(settings, vakint, cff_energy_degree_bound_reports)?;
        Ok(())
    }

    fn build_integrand(
        &mut self,
        settings: &GenerationSettings,
        vakint: &Vakint,
        cff_energy_degree_bound_reports: &Mutex<Vec<CffEnergyDegreeBoundReport>>,
    ) -> Result<TiVec<CutGroupId, ParametricIntegrands>> {
        let started = std::time::Instant::now();
        crate::debug_tags!(#generation, #profile, #uv, #graph, #summary;
            stage = "supergraph_build_integrand_start",
            graph = %self.graph.name,
            subtract_uv = settings.uv.subtract_uv,
            generate_integrated = settings.uv.generate_integrated,
            "Generation timing milestone"
        );
        let max_order = self
            .derived_data
            .cut_group_data
            .cut_groups
            .iter()
            .map(|cut_group| cut_group.related_esurface_group.max_occurence)
            .max()
            .unwrap();

        self.graph.param_builder.initialize_duals(max_order);
        crate::debug_tags!(#generation, #profile, #graph, #summary;
            stage = "supergraph_initialize_duals_done",
            graph = %self.graph.name,
            max_order,
            elapsed_ms = started.elapsed().as_secs_f64() * 1000.0,
            "Generation timing milestone"
        );

        let cuts = self
            .derived_data
            .cut_group_data
            .cut_groups
            .iter()
            .map(|cuts| CutSet {
                residue_selector: ResidueSelector {
                    lu: Some(cuts.lu_cut_selection(&self.graph, &self.cuts)),
                    left_th_cut: None,
                    right_th_cut: None,
                },
                union: cuts
                    .cuts
                    .iter()
                    .map(|cut_id| self.cuts[*cut_id].cut.as_subgraph())
                    .reduce(|cut_1, cut_2| cut_1.union(&cut_2))
                    .unwrap_or_else(|| self.graph.empty_subgraph()),
                canonicalize_external_shifts: false,
            })
            .collect();

        let cut_structure = CutStructure { cuts };
        crate::debug_tags!(#generation, #profile, #uv, #graph, #summary;
            stage = "supergraph_cutsets_done",
            graph = %self.graph.name,
            cut_count = cut_structure.cuts.len(),
            elapsed_ms = started.elapsed().as_secs_f64() * 1000.0,
            "Generation timing milestone"
        );

        let cff_options = self.graph.production_cff_3d_expression_options(settings)?;
        let production_expression = self
            .derived_data
            .global_cff_expression
            .as_ref()
            .expect("global_cff_expression should have been created");

        let orchestration_started = std::time::Instant::now();
        let parametric_integrands = settings.uv.orchestrator.parametric_integrands(
            &mut self.graph,
            cut_structure,
            vakint,
            OrientationProjection::exact_expression(
                production_expression,
                &cff_options,
                &settings.orientation_pattern,
                settings.explicit_orientation_sum_only,
            )
            .with_vacuum_subtraction(settings.medium.vacuum_subtraction)
            .with_energy_degree_bound_reports(cff_energy_degree_bound_reports),
            settings,
        )?;
        crate::debug_tags!(#generation, #profile, #uv, #graph, #summary;
            stage = "supergraph_parametric_orchestration_done",
            graph = %self.graph.name,
            parametric_integrand_count = parametric_integrands.len(),
            elapsed_ms = orchestration_started.elapsed().as_secs_f64() * 1000.0,
            total_elapsed_ms = started.elapsed().as_secs_f64() * 1000.0,
            "Generation timing milestone"
        );

        let finalize_started = std::time::Instant::now();
        let result = parametric_integrands
            .into_iter()
            .enumerate()
            .map(|(index, integrand)| {
                let group = &self.derived_data.cut_group_data.cut_groups[CutGroupId(index)];
                let lu_prefactor = group.lu_prefactor(&self.graph, &self.cuts)?;
                Ok(integrand.map(|a| a * &lu_prefactor))
            })
            .collect::<Result<_>>()?;
        crate::debug_tags!(#generation, #profile, #uv, #graph, #summary;
            stage = "supergraph_build_integrand_done",
            graph = %self.graph.name,
            elapsed_ms = finalize_started.elapsed().as_secs_f64() * 1000.0,
            total_elapsed_ms = started.elapsed().as_secs_f64() * 1000.0,
            "Generation timing milestone"
        );
        Ok(result)
    }

    fn threshold_rstar_power(subspace_loop_count: usize, is_on_right: bool) -> Atom {
        let radius_star = if is_on_right {
            GS.radius_star_right
        } else {
            GS.radius_star_left
        };
        Atom::var(radius_star).pow(subspace_loop_count as i32 * 3 - 1)
    }

    pub(crate) fn single_th_prefactor_helper_atoms(
        order: u8,
        subspace_loop_count: usize,
        is_on_right: bool,
        include_integrated: bool,
    ) -> (SingleThresholdPieces<Atom>, Atom) {
        let loop_3 = subspace_loop_count as i64 * 3;

        let i = Atom::i();

        let radius = if is_on_right {
            Atom::var(GS.radius_right)
        } else {
            Atom::var(GS.radius_left)
        };

        let radius_star = if is_on_right {
            Atom::var(GS.radius_star_right)
        } else {
            Atom::var(GS.radius_star_left)
        };
        let uv_damp_plus = if is_on_right {
            Atom::var(GS.uv_damp_plus_right)
        } else {
            Atom::var(GS.uv_damp_plus_left)
        };
        let uv_damp_minus = if is_on_right {
            Atom::var(GS.uv_damp_minus_right)
        } else {
            Atom::var(GS.uv_damp_minus_left)
        };
        let hfunction = if is_on_right {
            Atom::var(GS.hfunction_right_th)
        } else {
            Atom::var(GS.hfunction_left_th)
        };

        let laurent_coeff_indices = (1..=order).map(|i| -(i as i8));

        let mut laurent_coeffs = laurent_coeff_indices.map(|laurent_coeff_index| {
            build_derivative_structure_atom(order, laurent_coeff_index)
                .replace(GS.rescale_star)
                .with(radius_star.clone())
        });

        let delta_r_plus = &radius - &radius_star;
        let delta_r_minus = -&radius - &radius_star;

        let jacobian_ratio = (Atom::one() / &radius).pow(loop_3 - 1);

        let local_prefactor =
            &jacobian_ratio * (uv_damp_plus / &delta_r_plus + uv_damp_minus / &delta_r_minus);

        let integrated_prefactor = if include_integrated {
            // Runtime subtracts this helper. The left causal surface therefore
            // needs -i*pi here to restore +i*pi*delta(eta), and the right its conjugate.
            if is_on_right {
                i * Atom::var(GS.pi) * &jacobian_ratio * hfunction
            } else {
                -i * Atom::var(GS.pi) * &jacobian_ratio * hfunction
            }
        } else {
            Atom::zero()
        };

        let leading_laurent_coeff = laurent_coeffs.next().unwrap();
        // Keep the no-directive helper structurally identical to the historical single-output
        // lane. Evaluating the split pieces and adding them later changes f64 rounding.
        let mut legacy = (local_prefactor.clone() + integrated_prefactor.clone())
            * leading_laurent_coeff.clone();
        let mut local = local_prefactor * &leading_laurent_coeff;
        let integrated = integrated_prefactor * leading_laurent_coeff;

        for pow in 2..=order {
            let raised = laurent_coeffs.next().unwrap()
                * &jacobian_ratio
                * (Atom::one() / delta_r_plus.pow(pow as i64)
                    + Atom::one() / delta_r_minus.pow(pow as i64));
            legacy += raised.clone();
            local += raised;
        }

        debug!(
            "Threshold counterterm helper atoms for order {} and loop number {}: local={}, integrated={}, legacy={}",
            order, subspace_loop_count, local, integrated, legacy,
        );

        (SingleThresholdPieces { local, integrated }, legacy)
    }

    fn iterated_th_prefactor_helper_atoms(
        left_order: u8,
        right_order: u8,
        left_subspace_loop_count: usize,
        right_subspace_loop_count: usize,
        include_integrated: bool,
    ) -> (IteratedThresholdPieces<Atom>, Atom) {
        let (left_prefactors, left_legacy) = Self::single_th_prefactor_helper_atoms(
            left_order,
            left_subspace_loop_count,
            false,
            include_integrated,
        );
        let left_local = left_prefactors
            .local
            .replace(GS.eta)
            .with(Atom::var(GS.eta_left));
        let left_integrated = left_prefactors
            .integrated
            .replace(GS.eta)
            .with(Atom::var(GS.eta_left));

        let (right_prefactors, right_legacy) = Self::single_th_prefactor_helper_atoms(
            right_order,
            right_subspace_loop_count,
            true,
            include_integrated,
        );
        let right_local = right_prefactors
            .local
            .replace(GS.eta)
            .with(Atom::var(GS.eta_right));
        let right_integrated = right_prefactors
            .integrated
            .replace(GS.eta)
            .with(Atom::var(GS.eta_right));

        // The two single-side helpers start with the same generic `f` and `η` families.
        // Iterated subtraction evaluates one bivariate `f` and distinct left/right inverse-map
        // derivatives. Expansion exposes every product of univariate `f` derivatives so the
        // replacements below can fuse them into the corresponding bivariate derivative.
        let fuse = |left: Atom, right: Atom| {
            (left * right)
                .expand()
                .replace_multiple(Self::fuse_left_right_replacement())
        };

        let pieces = IteratedThresholdPieces {
            local_local: fuse(left_local.clone(), right_local.clone()),
            local_integrated: fuse(left_local, right_integrated.clone()),
            integrated_local: fuse(left_integrated.clone(), right_local),
            integrated_integrated: fuse(left_integrated, right_integrated),
        };
        let legacy = fuse(
            left_legacy.replace(GS.eta).with(Atom::var(GS.eta_left)),
            right_legacy.replace(GS.eta).with(Atom::var(GS.eta_right)),
        );

        (pieces, legacy)
    }

    pub(crate) fn single_th_prefactor_helper_params(order: u8, is_on_right: bool) -> Vec<Atom> {
        let radius_star = if is_on_right {
            Atom::var(GS.radius_star_right)
        } else {
            Atom::var(GS.radius_star_left)
        };

        let uv_damp_plus = if is_on_right {
            Atom::var(GS.uv_damp_plus_right)
        } else {
            Atom::var(GS.uv_damp_plus_left)
        };
        let uv_damp_minus = if is_on_right {
            Atom::var(GS.uv_damp_minus_right)
        } else {
            Atom::var(GS.uv_damp_minus_left)
        };
        let hfunction = if is_on_right {
            Atom::var(GS.hfunction_right_th)
        } else {
            Atom::var(GS.hfunction_left_th)
        };

        let radius = if is_on_right {
            Atom::var(GS.radius_right)
        } else {
            Atom::var(GS.radius_left)
        };

        let mut params = params_for_derivative_order(order)
            .into_iter()
            .map(|param| param.replace(GS.rescale_star).with(radius_star.clone()))
            .collect_vec();

        params.push(radius);
        params.push(radius_star);
        params.push(uv_damp_plus);
        params.push(uv_damp_minus);
        params.push(hfunction);
        params
    }

    pub(crate) fn iterated_th_prefactor_helper_params(
        left_order: u8,
        right_order: u8,
    ) -> Vec<Atom> {
        let mut iterated_params = params_for_iterated_threshold_ct(left_order, right_order);
        let left_radius_star = Atom::var(GS.radius_star_left);
        let right_radius_star = Atom::var(GS.radius_star_right);
        let left_radius = Atom::var(GS.radius_left);
        let right_radius = Atom::var(GS.radius_right);
        let left_uv_damp_plus = Atom::var(GS.uv_damp_plus_left);
        let left_uv_damp_minus = Atom::var(GS.uv_damp_minus_left);
        let right_uv_damp_plus = Atom::var(GS.uv_damp_plus_right);
        let right_uv_damp_minus = Atom::var(GS.uv_damp_minus_right);
        let left_hfunction = Atom::var(GS.hfunction_left_th);
        let right_hfunction = Atom::var(GS.hfunction_right_th);

        iterated_params.push(left_radius);
        iterated_params.push(left_radius_star);
        iterated_params.push(left_uv_damp_plus);
        iterated_params.push(left_uv_damp_minus);
        iterated_params.push(left_hfunction);
        iterated_params.push(right_radius);
        iterated_params.push(right_radius_star);
        iterated_params.push(right_uv_damp_plus);
        iterated_params.push(right_uv_damp_minus);
        iterated_params.push(right_hfunction);
        iterated_params
    }

    #[allow(clippy::too_many_arguments)]
    pub(crate) fn single_th_helper(
        &self,
        order: u8,
        subspace_loop_count: usize,
        is_on_right: bool,
        include_integrated: bool,
        outputs: LUThresholdHelperOutputs,
        dual_shape: Option<Vec<Vec<usize>>>,
        optimization_settings: OptimizationSettings,
        evaluator_settings: &EvaluatorSettings,
    ) -> Result<GenericEvaluator> {
        let (pieces, legacy) = Self::single_th_prefactor_helper_atoms(
            order,
            subspace_loop_count,
            is_on_right,
            include_integrated,
        );
        let atoms = match outputs {
            LUThresholdHelperOutputs::Legacy => vec![legacy],
            LUThresholdHelperOutputs::Pieces => vec![pieces.local, pieces.integrated],
            LUThresholdHelperOutputs::LegacyAndPieces => {
                vec![legacy, pieces.local, pieces.integrated]
            }
        };
        let params = Self::single_th_prefactor_helper_params(order, is_on_right);

        let fn_map = FunctionMap::new();

        let evaluator = GenericEvaluator::new_from_raw_params(
            atoms,
            &params,
            &fn_map,
            vec![],
            optimization_settings,
            dual_shape.map(|shape| (shape, Vec::new())),
            evaluator_settings,
        )?
        .into_eager_only();

        Ok(evaluator)
    }

    #[allow(clippy::too_many_arguments)]
    pub(crate) fn iterated_th_helper(
        &self,
        left_order: u8,
        right_order: u8,
        left_subspace_loop_count: usize,
        right_subspace_loop_count: usize,
        include_integrated: bool,
        outputs: LUThresholdHelperOutputs,
        dual_shape: Option<Vec<Vec<usize>>>,
        optimization_settings: OptimizationSettings,
        evaluator_settings: &EvaluatorSettings,
    ) -> Result<GenericEvaluator> {
        let (pieces, legacy) = Self::iterated_th_prefactor_helper_atoms(
            left_order,
            right_order,
            left_subspace_loop_count,
            right_subspace_loop_count,
            include_integrated,
        );
        let atoms = match outputs {
            LUThresholdHelperOutputs::Legacy => vec![legacy],
            LUThresholdHelperOutputs::Pieces => vec![
                pieces.local_local,
                pieces.local_integrated,
                pieces.integrated_local,
                pieces.integrated_integrated,
            ],
            LUThresholdHelperOutputs::LegacyAndPieces => vec![
                legacy,
                pieces.local_local,
                pieces.local_integrated,
                pieces.integrated_local,
                pieces.integrated_integrated,
            ],
        };

        let params = Self::iterated_th_prefactor_helper_params(left_order, right_order);

        let fn_map = FunctionMap::new();

        let evaluator = GenericEvaluator::new_from_raw_params(
            atoms,
            &params,
            &fn_map,
            vec![],
            optimization_settings,
            dual_shape.map(|shape| (shape, Vec::new())),
            evaluator_settings,
        )?
        .into_eager_only();

        Ok(evaluator)
    }

    fn fuse_left_right_replacement() -> Vec<Replacement> {
        let f = symbol!("f");

        vec![
            Replacement::new(
                (function!(f, GS.radius_star_left) * function!(f, GS.radius_star_right))
                    .to_pattern(),
                function!(f, GS.radius_star_left, GS.radius_star_right),
            ),
            Replacement::new(
                (function!(f, GS.radius_star_left)
                    * function!(Symbol::DERIVATIVE, W_.x_, f, GS.radius_star_right))
                .to_pattern(),
                function!(
                    Symbol::DERIVATIVE,
                    0,
                    W_.x_,
                    f,
                    GS.radius_star_left,
                    GS.radius_star_right
                ),
            ),
            Replacement::new(
                (function!(Symbol::DERIVATIVE, W_.x_, f, GS.radius_star_left)
                    * function!(f, GS.radius_star_right))
                .to_pattern(),
                function!(
                    Symbol::DERIVATIVE,
                    W_.x_,
                    0,
                    f,
                    GS.radius_star_left,
                    GS.radius_star_right
                ),
            ),
            Replacement::new(
                (function!(Symbol::DERIVATIVE, W_.x_, f, GS.radius_star_left)
                    * function!(Symbol::DERIVATIVE, W_.y_, f, GS.radius_star_right))
                .to_pattern(),
                function!(
                    Symbol::DERIVATIVE,
                    W_.x_,
                    W_.y_,
                    f,
                    GS.radius_star_left,
                    GS.radius_star_right
                ),
            ),
        ]
    }
    //fn th_prefactor_helper(
    //    &self,
    //    subspace_loop_count: usize,
    //    is_on_right: bool,
    //    include_integrated: bool,
    //) -> Atom {
    //    let radius = if is_on_right {
    //        Atom::var(GS.radius_right)
    //    } else {
    //        Atom::var(GS.radius_left)
    //    };

    //    let radius_star = if is_on_right {
    //        Atom::var(GS.radius_star_right)
    //    } else {
    //        Atom::var(GS.radius_star_left)
    //    };

    //    let grad_eta = if is_on_right {
    //        Atom::var(GS.deta_right_th)
    //    } else {
    //        Atom::var(GS.deta_left_th)
    //    };

    //    let uv_damp_plus = if is_on_right {
    //        Atom::var(GS.uv_damp_plus_right)
    //    } else {
    //        Atom::var(GS.uv_damp_plus_left)
    //    };

    //    let uv_damp_minus = if is_on_right {
    //        Atom::var(GS.uv_damp_minus_right)
    //    } else {
    //        Atom::var(GS.uv_damp_minus_left)
    //    };

    //    let h_function = if is_on_right {
    //        Atom::var(GS.hfunction_right_th)
    //    } else {
    //        Atom::var(GS.hfunction_left_th)
    //    };

    //    let i = if is_on_right { Atom::i() } else { -Atom::i() };
    //    let pi = Atom::var(GS.pi);

    //    let jacobian_ratio = (&radius_star / &radius).pow(subspace_loop_count as i64 * 3 - 1);

    //    let local_prefactor = (&uv_damp_plus / (&radius - &radius_star)
    //        + &uv_damp_minus / (-&radius - &radius_star))
    //        / &grad_eta
    //        * &jacobian_ratio;

    //    let integrated_prefactor = if include_integrated {
    //        &h_function * &i * &pi * &jacobian_ratio / &grad_eta
    //    } else {
    //        Atom::new()
    //    };

    //    debug!(
    //        "th prefactor local: {}, integrated: {}",
    //        local_prefactor, integrated_prefactor
    //    );

    //    local_prefactor + integrated_prefactor
    //}

    fn build_lmbs(&mut self) -> Result<()> {
        let mut lmbs: TiVec<LmbIndex, LoopMomentumBasis> = vec![].into();

        let externals: SuBitGraph = self.graph.empty_subgraph();
        let full_filter = self.graph.full_filter();
        let cut_graph = full_filter.subtract(&self.graph.initial_state_cut.right);

        for s in self.graph.all_spanning_forests_of(&cut_graph) {
            let mut lmb = self.graph.lmb_impl(&full_filter, &s, externals.clone())?;
            let mut exts = vec![];

            for i in lmb.loop_edges.iter() {
                let (_, p) = &self.graph[i];

                if self.graph.initial_state_cut.intersects(p) {
                    exts.push(*i);
                }
            }

            exts.sort();

            for e in exts {
                let mut loopid = None;
                for (l, s) in lmb.edge_signatures[e].internal.iter_enumerated() {
                    if s.is_non_zero() {
                        if loopid.is_none() {
                            loopid = Some(l);
                        } else {
                            panic!("external edge has multiple loop momenta")
                        }
                    }
                }
                lmb.put_loop_to_ext(loopid.unwrap());
            }
            let external_momentum_edge_order = self.graph.external_momentum_edge_order();
            lmb.canonicalize_external_order(&external_momentum_edge_order);
            //lmbs.push(self.graph.lmb_impl(&full_filter, &s, externals.clone()));
            lmbs.push(lmb);
        }

        let sorted_graph_loop_edges = self
            .graph
            .loop_momentum_basis
            .loop_edges
            .iter()
            .copied()
            .sorted()
            .collect_vec();

        let matching_lmb_index = lmbs.iter_enumerated().find_map(|(lmb_index, lmb)| {
            let sorted_loop_edges = lmb.loop_edges.iter().copied().sorted().collect_vec();
            (sorted_loop_edges == sorted_graph_loop_edges).then_some(lmb_index)
        });

        if let Some(matching_lmb_index) = matching_lmb_index {
            {
                let matching_lmb = &mut lmbs[matching_lmb_index];
                for (target_pos, target_edge) in
                    self.graph.loop_momentum_basis.loop_edges.iter_enumerated()
                {
                    let current_pos = matching_lmb
                        .loop_edges
                        .iter()
                        .find_position(|edge| *edge == target_edge)
                        .map(|(pos, _)| pos)
                        .unwrap();

                    if current_pos != target_pos.0 {
                        matching_lmb.swap_loops(current_pos.into(), target_pos);
                    }
                }
            }

            let matching_lmb_pos: usize = matching_lmb_index.into();
            if matching_lmb_pos != 0 {
                lmbs.raw[..=matching_lmb_pos].rotate_right(1);
            }
        } else {
            warn!(
                "Could not match current graph LMB against generated LMBs for graph {}",
                self.graph.name
            );
        }

        let sorted_graph_loop_edges = self
            .graph
            .loop_momentum_basis
            .loop_edges
            .iter()
            .copied()
            .sorted()
            .collect_vec();

        let matching_lmb_index = lmbs.iter_enumerated().find_map(|(lmb_index, lmb)| {
            let sorted_loop_edges = lmb.loop_edges.iter().copied().sorted().collect_vec();
            (sorted_loop_edges == sorted_graph_loop_edges).then_some(lmb_index)
        });

        if let Some(matching_lmb_index) = matching_lmb_index {
            {
                let matching_lmb = &mut lmbs[matching_lmb_index];
                for (target_pos, target_edge) in
                    self.graph.loop_momentum_basis.loop_edges.iter_enumerated()
                {
                    let current_pos = matching_lmb
                        .loop_edges
                        .iter()
                        .find_position(|edge| *edge == target_edge)
                        .map(|(pos, _)| pos)
                        .unwrap();

                    if current_pos != target_pos.0 {
                        matching_lmb.swap_loops(current_pos.into(), target_pos);
                    }
                }
            }

            let matching_lmb_pos: usize = matching_lmb_index.into();
            if matching_lmb_pos != 0 {
                lmbs.raw[..=matching_lmb_pos].rotate_right(1);
            }
        } else {
            warn!(
                "Could not match current graph LMB against generated LMBs for graph {}",
                self.graph.name
            );
        }

        self.derived_data.lmbs = Some(lmbs);
        Ok(())
    }

    fn build_multi_channeling_channels(&mut self, override_lmb_heuristics: bool) -> Result<()> {
        let lmbs = self.derived_data.lmbs.as_ref().unwrap();
        let channels = if override_lmb_heuristics {
            self.graph
                .build_multi_channeling_channels(lmbs, override_lmb_heuristics)
        } else {
            self.build_cross_section_multi_channeling_channels(lmbs)?
        };

        self.derived_data.multi_channeling_setup = Some(channels);
        Ok(())
    }

    fn build_cross_section_multi_channeling_channels(
        &self,
        lmbs: &TiVec<LmbIndex, LoopMomentumBasis>,
    ) -> Result<LmbMultiChannelingSetup> {
        let channels = self.select_cross_section_lmb_channel_indices(lmbs)?;
        debug!(
            "number of lmbs: {}, number of cross-section channels: {}",
            lmbs.len(),
            channels.len()
        );
        Ok(LmbMultiChannelingSetup {
            master_edge_masses: Default::default(),
            sampling_bridge: Default::default(),
            sampling_bridge_quad: Default::default(),
            sampling_bridge_fixed256: Default::default(),
            sampling_bridge_arb: Default::default(),
            sampling_source: Default::default(),
            sampling_catalogue: Default::default(),
            sampling_programs: Default::default(),
            lmb_basis_ids: channels,
            graph: self.graph.clone(),
            all_bases: lmbs.clone(),
        })
    }

    fn select_cross_section_lmb_channel_indices(
        &self,
        lmbs: &TiVec<LmbIndex, LoopMomentumBasis>,
    ) -> Result<TiVec<crate::integrands::process::SamplingChannelId, LmbIndex>> {
        if self.cuts.is_empty() {
            return Ok(self.graph.select_amplitude_lmb_channel_indices(
                lmbs,
                false,
                LmbChannelFallback::CurrentGraphBasis,
            ));
        }

        let mut lmb_index_by_loop_edges = BTreeMap::<Vec<EdgeIndex>, LmbIndex>::new();
        for (lmb_index, lmb) in lmbs.iter_enumerated() {
            lmb_index_by_loop_edges
                .entry(lmb.loop_edges.iter().copied().sorted().collect())
                .or_insert(lmb_index);
        }

        let mut selected_loop_edge_sets = BTreeSet::<Vec<EdgeIndex>>::new();
        for cut in &self.cuts {
            let cut_edges = cut.cut.as_subgraph();
            let cut_edge_ids = self
                .graph
                .underlying
                .iter_edges_of(&cut_edges)
                .map(|(_, edge_id, _)| edge_id)
                .sorted()
                .collect_vec();
            let Some(excluded_cut_edge) = self.excluded_cut_edge_for_lmb_channel(&cut_edge_ids)
            else {
                continue;
            };

            let fixed_cut_edges = cut_edge_ids
                .iter()
                .copied()
                .filter(|edge_id| *edge_id != excluded_cut_edge)
                .collect::<BTreeSet<_>>();
            let left_basis_edge_sets =
                self.selected_cut_side_lmb_edge_sets(cut.left.subtract(&cut_edges))?;
            let right_basis_edge_sets =
                self.selected_cut_side_lmb_edge_sets(cut.right.subtract(&cut_edges))?;

            for left_basis_edges in &left_basis_edge_sets {
                for right_basis_edges in &right_basis_edge_sets {
                    let loop_edges = fixed_cut_edges
                        .iter()
                        .copied()
                        .chain(left_basis_edges.iter().copied())
                        .chain(right_basis_edges.iter().copied())
                        .collect::<BTreeSet<_>>()
                        .into_iter()
                        .collect_vec();
                    selected_loop_edge_sets.insert(loop_edges);
                }
            }
        }

        let mut channels = Vec::<LmbIndex>::new();
        for loop_edges in selected_loop_edge_sets {
            let lmb_index = lmb_index_by_loop_edges.get(&loop_edges).ok_or_else(|| {
                eyre!(
                    "Could not find a generated cross-section LMB with loop edges [{}] for graph '{}'.",
                    loop_edges.iter().map(|edge| edge.to_string()).join(", "),
                    self.graph.name
                )
            })?;
            channels.push(*lmb_index);
        }

        if channels.is_empty() {
            Ok(self.graph.select_amplitude_lmb_channel_indices(
                lmbs,
                false,
                LmbChannelFallback::CurrentGraphBasis,
            ))
        } else {
            Ok(channels.into_iter().sorted().dedup().collect())
        }
    }

    fn excluded_cut_edge_for_lmb_channel(&self, cut_edge_ids: &[EdgeIndex]) -> Option<EdgeIndex> {
        Self::excluded_cut_edge_for_lmb_channel_in(&self.graph, cut_edge_ids)
    }

    fn excluded_cut_edge_for_lmb_channel_in(
        graph: &Graph,
        cut_edge_ids: &[EdgeIndex],
    ) -> Option<EdgeIndex> {
        cut_edge_ids
            .iter()
            .copied()
            .filter(|edge_id| graph[*edge_id].particle.is_massive())
            .min()
            .or_else(|| {
                cut_edge_ids
                    .iter()
                    .copied()
                    .filter(|edge_id| graph[*edge_id].particle.is_fermion())
                    .min()
            })
            .or_else(|| cut_edge_ids.iter().copied().min())
    }

    fn selected_cut_side_lmb_edge_sets(
        &self,
        mut side_subgraph: SuBitGraph,
    ) -> Result<Vec<Vec<EdgeIndex>>> {
        for (pair, _, edge) in self.graph.underlying.iter_edges() {
            if edge.data.is_dummy {
                side_subgraph.sub(pair);
            }
        }

        if self.graph.underlying.cyclotomatic_number(&side_subgraph) == 0 {
            return Ok(vec![Vec::new()]);
        }

        let side_lmbs = self.graph.generate_loop_momentum_bases_of(&side_subgraph);
        if side_lmbs.is_empty() {
            return Err(eyre!(
                "Could not generate cut-side LMBs for a non-tree side of graph '{}'.",
                self.graph.name
            ));
        }

        let selected_side_lmbs = self.graph.select_amplitude_lmb_channel_indices(
            &side_lmbs,
            false,
            LmbChannelFallback::FirstBasis,
        );
        let mut edge_sets = BTreeSet::<Vec<EdgeIndex>>::new();
        for lmb_index in selected_side_lmbs {
            edge_sets.insert(
                side_lmbs[lmb_index]
                    .loop_edges
                    .iter()
                    .copied()
                    .sorted()
                    .collect(),
            );
        }

        if edge_sets.is_empty() {
            Ok(vec![Vec::new()])
        } else {
            Ok(edge_sets.into_iter().collect())
        }
    }

    fn threshold_counterterm_association(
        &self,
        threshold_candidate: &TopologicalThresholdCandidate,
        sandwich: &SuBitGraph,
        cut: &OrientedCut,
    ) -> ThresholdCountertermAssociation {
        let boundary_edges = |boundary: &OrientedCut| {
            let boundary_filter = boundary.left.union(&boundary.right).intersection(sandwich);
            let edges = self
                .graph
                .iter_edges_of(&boundary_filter)
                .map(|(_, edge_id, _)| edge_id)
                .sorted()
                .collect();
            let has_left_orientation =
                boundary_filter.intersection(&boundary.left).n_included() > 0;
            let has_right_orientation =
                boundary_filter.intersection(&boundary.right).n_included() > 0;
            (edges, has_left_orientation != has_right_orientation)
        };

        let (cut_boundary_edges, cut_bound_is_applicable) = boundary_edges(cut);
        let (threshold_boundary_edges, threshold_bound_is_applicable) =
            boundary_edges(&threshold_candidate.cut);

        ThresholdCountertermAssociation {
            topological_threshold_id: threshold_candidate.topological_threshold_id,
            esurface_id: threshold_candidate
                .selected_orientation_esurface_id
                .expect("only selected-orientation threshold instances have associations"),

            cut_boundary_edges,
            threshold_boundary_edges,
            invariant_bound_is_applicable: cut_bound_is_applicable && threshold_bound_is_applicable,
        }
    }

    #[allow(clippy::too_many_arguments)]
    fn classify_threshold_counterterm_association(
        &self,
        association: &ThresholdCountertermAssociation,
        cut_id: CutId,
        cut: &CrossSectionCut,
        subspace: &SubspaceData,
        all_lmbs: &TiVec<LmbIndex, LoopMomentumBasis>,
        model: &Model,
        runtime_settings: &RuntimeSettings,
        settings: &GenerationSettings,
    ) -> ThresholdCountertermStatus {
        let esurface = &self.derived_data.topological_threshold_esurfaces
            [association.topological_threshold_id];
        if !esurface.has_radial_dependence_in_subspace(subspace, all_lmbs, &self.graph) {
            debug!(
                "Skipping graph '{}' cut {} threshold E-surface {} after cut-relative classification {:?}",
                self.graph.name,
                cut_id.0,
                association.topological_threshold_id.0,
                ThresholdCountertermStatus::NoRadialDependence,
            );
            return ThresholdCountertermStatus::NoRadialDependence;
        }

        let structural_status = if association.invariant_bound_is_applicable {
            self.graph.classify_threshold_pinch(
                &association.cut_boundary_edges,
                &association.threshold_boundary_edges,
            )
        } else {
            ThresholdPinchStatus::NotProven
        };
        if structural_status == ThresholdPinchStatus::Always {
            debug!(
                "Skipping graph '{}' cut {} threshold E-surface {} after cut-relative classification {:?}",
                self.graph.name,
                cut_id.0,
                association.topological_threshold_id.0,
                ThresholdCountertermStatus::AlwaysPinched,
            );
            return ThresholdCountertermStatus::AlwaysPinched;
        }

        let status = if settings.threshold_subtraction.check_esurface_at_generation {
            association.classify_for_model(
                &self.graph,
                cut,
                model,
                &self.graph.param_builder,
                runtime_settings,
                settings.threshold_subtraction.esurface_existence_threshold,
            )
        } else {
            match structural_status {
                ThresholdPinchStatus::CanBecome => ThresholdCountertermStatus::CanBecomePinched,
                ThresholdPinchStatus::NotProven => ThresholdCountertermStatus::PotentiallyExisting,
                ThresholdPinchStatus::Always => unreachable!(
                    "identically pinched threshold associations return before classification"
                ),
            }
        };
        if !status
            .is_eligible_for_generation(settings.threshold_subtraction.check_esurface_at_generation)
        {
            debug!(
                "Skipping graph '{}' cut {} threshold E-surface {} after cut-relative classification {:?}",
                self.graph.name, cut_id.0, association.topological_threshold_id.0, status,
            );
        }
        status
    }

    fn topological_threshold_candidates(
        &mut self,
        settings: &GenerationSettings,
    ) -> Result<Vec<TopologicalThresholdCandidate>> {
        let global_cff = self
            .derived_data
            .global_cff_expression
            .as_ref()
            .expect("global_cff_expression should have been created");
        // CFF generation retains the full catalogue for 3D UV construction.
        // Only requested orientations can activate threshold-counterterm instances.
        let selected_orientations = global_cff
            .expression
            .orientations
            .iter()
            .filter(|orientation| settings.orientation_pattern.filter(*orientation))
            .collect_vec();
        let selected_cff_esurface_ids = selected_orientations
            .iter()
            .flat_map(|orientation| {
                orientation
                    .iter_denominator_nodes()
                    .filter_map(|node| match node.data {
                        HybridSurfaceID::Esurface(esurface_id) => Some(esurface_id),
                        _ => None,
                    })
            })
            .collect::<BTreeSet<_>>();

        let mut candidates = self.graph.all_st_cuts_for_cs(
            self.source_nodes.clone(),
            self.target_nodes.clone(),
            &self.graph.get_initial_state_tree(),
        );
        candidates.retain(|(_left, cut, _right)| cut.nedges(&self.graph) > 1);
        candidates.sort_by(|a, b| a.1.cmp(&b.1));

        let mut topological_threshold_esurfaces = TiVec::new();
        let candidates: Vec<Result<TopologicalThresholdCandidate>> = candidates
            .into_iter()
            .map(|(left, cut, right)| {
                let threshold_esurface = Esurface::new_from_cut_left(
                    &self.graph.underlying,
                    &CrossSectionCut {
                        cut: cut.clone(),
                        left: left.clone(),
                        right: right.clone(),
                    },
                    Some(&self.graph.initial_state_cut),
                );
                let has_selected_orientation_instance =
                    selected_orientations.iter().any(|orientation| {
                        let selected_orientation = &orientation.data.orientation;
                        self.graph
                            .iter_edges_of(&cut)
                            .map(|(_, edge_id, _)| edge_id)
                            .zip(cut.iter_edges(&self.graph).map(|(direction, _)| direction))
                            .all(|(edge_id, direction)| {
                                selected_orientation[edge_id] == direction
                            })
                    });
                if !has_selected_orientation_instance {
                    debug!(
                        "Keeping topology-discovered threshold surface {:?} dormant in graph '{}': no selected-orientation instance occurs",
                        threshold_esurface.energies, self.graph.name,
                    );
                    let topological_threshold_id = TopologicalThresholdId(
                        topological_threshold_esurfaces.len(),
                    );
                    topological_threshold_esurfaces.push(threshold_esurface.clone());
                    return Ok(TopologicalThresholdCandidate {
                        topological_threshold_id,
                        left,
                        cut,
                        right,
                        threshold_esurface,
                        selected_orientation_esurface_id: None,
                    });
                }

                let esurface_id = match self
                    .graph
                    .surface_cache
                    .esurface_cache
                    .position(|esurface| esurface == &threshold_esurface)
                {
                    Some(esurface_id) => esurface_id,
                    None => {
                        return Err(eyre!(
                            "Topology-discovered threshold surface {:?} has a selected-orientation instance but is missing from graph '{}' selected CFF surface cache",
                            threshold_esurface.energies,
                            self.graph.name,
                        ));
                    }
                };
                if !selected_cff_esurface_ids.contains(&esurface_id) {
                    return Err(eyre!(
                        "Topology-discovered threshold surface {:?} has a selected-orientation instance in graph '{}' but E-surface {} does not occur in its selected CFF expression",
                        threshold_esurface.energies,
                        self.graph.name,
                        esurface_id.0,
                    ));
                }

                let topological_threshold_id = TopologicalThresholdId(
                    topological_threshold_esurfaces.len(),
                );
                topological_threshold_esurfaces.push(threshold_esurface.clone());
                Ok(TopologicalThresholdCandidate {
                    topological_threshold_id,
                    left,
                    cut,
                    right,
                    threshold_esurface,
                    selected_orientation_esurface_id: Some(esurface_id),
                })

            })
            .collect();
        let candidates = candidates.into_iter().collect::<Result<Vec<_>>>()?;
        self.derived_data.topological_threshold_esurfaces = topological_threshold_esurfaces;
        Ok(candidates)
    }

    fn physical_cut_edges(&self, cut_id: CutId) -> Vec<EdgeIndex> {
        self.graph
            .iter_edges_of(&self.cuts[cut_id].cut)
            .map(|(_, edge_id, _)| edge_id)
            .sorted()
            .collect()
    }

    #[cfg(test)]
    fn threshold_energy_edges(&self, esurface_id: EsurfaceID) -> Vec<EdgeIndex> {
        self.graph.surface_cache.esurface_cache[esurface_id]
            .energies
            .iter()
            .copied()
            .sorted()
            .collect()
    }

    fn cut_group_side_subgraphs(
        &self,
        cut_group_id: CutGroupId,
    ) -> Result<(SuBitGraph, SuBitGraph)> {
        let cut_group = &self.derived_data.cut_group_data.cut_groups[cut_group_id];
        let smallest_side = |side: ThresholdCountertermSide| {
            cut_group
                .cuts
                .iter()
                .map(|cut_id| match side {
                    ThresholdCountertermSide::Amplitude => {
                        unreachable!("amplitude side is not valid in cross-section resolution")
                    }
                    ThresholdCountertermSide::Left => self.cuts[*cut_id].left.clone(),
                    ThresholdCountertermSide::Right => self.cuts[*cut_id].right.clone(),
                })
                .min_by_key(|subgraph| subgraph.nedges(&self.graph))
                .ok_or_else(|| {
                    eyre!(
                        "Graph '{}' cut group {} has no {:?} subgraph",
                        self.graph.name,
                        cut_group_id.0,
                        side,
                    )
                })
        };

        Ok((
            smallest_side(ThresholdCountertermSide::Left)?,
            smallest_side(ThresholdCountertermSide::Right)?,
        ))
    }

    fn cut_group_compatible_parent_lmbs(&self, cut_group_id: CutGroupId) -> Result<Vec<LmbIndex>> {
        let representative_cut_id = self.derived_data.cut_group_data.cut_groups[cut_group_id]
            .cuts
            .first()
            .copied()
            .ok_or_else(|| {
                eyre!(
                    "Graph '{}' has an empty cut group {} while selecting threshold-counterterm parents",
                    self.graph.name,
                    cut_group_id.0,
                )
            })?;
        let all_lmbs = self
            .derived_data
            .lmbs
            .as_ref()
            .expect("threshold resolution requires loop-momentum bases");
        Ok(all_lmbs
            .iter_enumerated()
            .filter_map(|(lmb_index, lmb)| {
                let mut cut_edges = self
                    .graph
                    .iter_edges_of(&self.cuts[representative_cut_id].cut)
                    .map(|(_, edge_id, _)| edge_id)
                    .collect_vec();
                cut_edges.retain(|edge_id| !lmb.loop_edges.contains(edge_id));
                (cut_edges.len() == 1).then_some(lmb_index)
            })
            .collect())
    }

    fn configured_threshold_variants(
        spec: &crate::graph::threshold_counterterms::ThresholdCountertermSpec,
        cut_edges: &[EdgeIndex],
        threshold_edges: &[EdgeIndex],
    ) -> Vec<(ThresholdCountertermVariant, ThresholdCountertermOrigin)> {
        let configured = spec
            .cuts
            .iter()
            .find(|cut| cut.edges == cut_edges)
            .and_then(|cut| {
                cut.thresholds
                    .iter()
                    .find(|threshold| threshold.edges == threshold_edges)
            })
            .filter(|threshold| !threshold.counterterms.is_empty());

        match configured {
            Some(threshold) => threshold
                .counterterms
                .iter()
                .cloned()
                .map(|mut variant| {
                    if let Some(multiplier) = &mut variant.multiplier {
                        for (name, definition) in &spec.function_map {
                            multiplier
                                .function_map
                                .entry(name.clone())
                                .or_insert_with(|| definition.clone());
                        }
                    }
                    (variant, ThresholdCountertermOrigin::Explicit)
                })
                .collect(),
            None => vec![(
                ThresholdCountertermVariant {
                    name: Some("default".to_string()),
                    subspace: None,
                    parent_lmb: None,
                    disable: false,
                    multiplier: None,
                    group_id: None,
                },
                ThresholdCountertermOrigin::Autogenerated,
            )],
        }
    }

    fn resolve_threshold_variant_subspace(
        &self,
        variant: &ThresholdCountertermVariant,
        legacy_subspace: &SubspaceData,
        containing_subgraph: &SuBitGraph,
        context: &str,
    ) -> Result<SubspaceData> {
        let Some(requested_parent_lmb) = &variant.parent_lmb else {
            if variant.subspace.is_some() {
                return Err(eyre!(
                    "{context} requires a full parent_lmb for its explicit subspace"
                ));
            }
            return Ok(legacy_subspace.clone());
        };
        let all_lmbs = self
            .derived_data
            .lmbs
            .as_ref()
            .expect("threshold resolution requires loop-momentum bases");
        let parent_lmb_index = all_lmbs
            .iter_enumerated()
            .find_map(|(lmb_index, lmb)| {
                lmb.loop_edges
                    .iter()
                    .eq(requested_parent_lmb.iter())
                    .then_some(lmb_index)
            })
            .ok_or_else(|| {
                eyre!(
                    "{context} requests parent_lmb {:?}, which is not among graph '{}' generated full LMBs",
                    requested_parent_lmb,
                    self.graph.name,
                )
            })?;

        match &variant.subspace {
            Some(requested_basis_edges) => SubspaceData::new_from_parent_basis_edges(
                requested_basis_edges,
                containing_subgraph,
                parent_lmb_index,
                &self.graph,
                all_lmbs,
            )
            .or_else(|side_error| {
                SubspaceData::new_from_parent_basis_edges(
                    requested_basis_edges,
                    &self.graph.full_filter(),
                    parent_lmb_index,
                    &self.graph,
                    all_lmbs,
                )
                .with_context(|| {
                    format!("side-projected embedding was rejected with: {side_error:#}")
                })
            }),
            None => SubspaceData::new_with_user_selected_lmb(
                containing_subgraph.clone(),
                parent_lmb_index,
                &self.graph,
                all_lmbs,
            ),
        }
        .with_context(|| format!("{context} is incompatible with its requested parent_lmb"))
    }

    fn threshold_subspaces_are_equivalent(
        left: &SubspaceData,
        right: &SubspaceData,
        all_lmbs: &TiVec<LmbIndex, LoopMomentumBasis>,
    ) -> bool {
        left.has_equivalent_embedding(right, all_lmbs)
    }

    fn resolve_threshold_counterterm_directives(
        &self,
        model: &Model,
        settings: &GenerationSettings,
        runtime_settings: &RuntimeSettings,
        all_possible_thresholds: &[TopologicalThresholdCandidate],
    ) -> Result<(
        ResolvedThresholdCounterterms,
        TiVec<CutId, CutThresholdCountertermAssociations>,
    )> {
        let all_lmbs = self
            .derived_data
            .lmbs
            .as_ref()
            .expect("threshold resolution requires loop-momentum bases");
        let mut cut_group_by_cut = vec![None; self.cuts.len()];
        for (cut_group_id, cut_group) in self
            .derived_data
            .cut_group_data
            .cut_groups
            .iter_enumerated()
        {
            for cut_id in &cut_group.cuts {
                cut_group_by_cut[cut_id.0] = Some(cut_group_id);
            }
        }

        let threshold_raised_data = self.graph.determine_raised_esurfaces_from_expression(
            &self
                .derived_data
                .global_cff_expression
                .as_ref()
                .expect("global_cff_expression should have been created")
                .expression,
        );
        let mut raised_threshold_ids: TiVec<EsurfaceID, Option<RaisedEsurfaceId>> =
            ti_vec![None; self.graph.surface_cache.esurface_cache.len()];
        for (raised_threshold_id, raised_group) in
            threshold_raised_data.raised_groups.iter_enumerated()
        {
            for &esurface_id in &raised_group.esurface_ids {
                raised_threshold_ids[esurface_id] = Some(raised_threshold_id);
            }
        }

        let mut physical_drafts = Vec::new();
        let mut cut_threshold_associations: TiVec<CutId, CutThresholdCountertermAssociations> =
            ti_vec![CutThresholdCountertermAssociations::default(); self.cuts.len()];
        let mut legacy_equivalent = true;
        let mut matched_explicit_cuts = BTreeSet::new();
        let mut matched_explicit_thresholds = BTreeSet::new();
        let mut explicit_group_members = BTreeMap::<usize, Vec<_>>::new();

        for (cut_id, cut) in self.cuts.iter_enumerated() {
            let cut_group_id = cut_group_by_cut[cut_id.0].ok_or_else(|| {
                eyre!(
                    "Physical cut {} in graph '{}' is missing from the cut-group partition",
                    cut_id.0,
                    self.graph.name,
                )
            })?;
            let cut_edges = self.physical_cut_edges(cut_id);
            if self
                .graph
                .threshold_counterterms
                .cuts
                .iter()
                .any(|configured| configured.edges == cut_edges)
            {
                matched_explicit_cuts.insert(cut_edges.clone());
            }
            let (left_containing_subgraph, right_containing_subgraph) =
                self.cut_group_side_subgraphs(cut_group_id)?;
            let (legacy_left_subspace, legacy_right_subspace) =
                &self.derived_data.subspace_data[cut_group_id];

            for threshold_candidate in all_possible_thresholds {
                let inferred = if cut.left.includes(&threshold_candidate.left) {
                    let sandwich = cut.left.intersection(&threshold_candidate.right);
                    self.graph.underlying.is_connected(&sandwich).then_some((
                        ThresholdCountertermSide::Left,
                        sandwich,
                        legacy_left_subspace,
                        &left_containing_subgraph,
                    ))
                } else if cut.right.includes(&threshold_candidate.right) {
                    let sandwich = cut.right.intersection(&threshold_candidate.left);
                    self.graph.underlying.is_connected(&sandwich).then_some((
                        ThresholdCountertermSide::Right,
                        sandwich,
                        legacy_right_subspace,
                        &right_containing_subgraph,
                    ))
                } else {
                    None
                };
                let Some((side, sandwich, legacy_subspace, containing_subgraph)) = inferred else {
                    continue;
                };

                let threshold_edges = threshold_candidate
                    .threshold_esurface
                    .energies
                    .iter()
                    .copied()
                    .sorted()
                    .collect_vec();
                if self
                    .graph
                    .threshold_counterterms
                    .cuts
                    .iter()
                    .find(|configured| configured.edges == cut_edges)
                    .is_some_and(|configured| {
                        configured
                            .thresholds
                            .iter()
                            .any(|threshold| threshold.edges == threshold_edges)
                    })
                {
                    matched_explicit_thresholds
                        .insert((cut_edges.clone(), threshold_edges.clone()));
                }
                // Explicit groups also constrain dormant, disabled and generation-ineligible
                // associations. Resolve their declared cycles before orientation/runtime filters.
                for (variant, _) in Self::configured_threshold_variants(
                    &self.graph.threshold_counterterms,
                    &cut_edges,
                    &threshold_edges,
                ) {
                    let Some(group_id) = variant.group_id else {
                        continue;
                    };
                    let context = format!(
                        "graph '{}' cut {:?} threshold {:?} variant '{}'",
                        self.graph.name,
                        cut_edges.iter().map(|edge| edge.0).collect_vec(),
                        threshold_edges.iter().map(|edge| edge.0).collect_vec(),
                        variant.name.as_deref().unwrap_or("default"),
                    );
                    let subspace = self.resolve_threshold_variant_subspace(
                        &variant,
                        legacy_subspace,
                        containing_subgraph,
                        &context,
                    )?;
                    let signature = subspace.solve_signature(all_lmbs);
                    let signed_cycles = signature
                        .iter()
                        .map(|(edge, cycle)| {
                            (
                                edge.0,
                                cycle
                                    .iter()
                                    .map(|(edge, sign)| (edge.0, sign))
                                    .collect_vec(),
                            )
                        })
                        .collect_vec();
                    let description = format!(
                        "  {context}, cut_group={}, cut={}, side={side:?}, parent={:?}, signed_cycles={signed_cycles:?}, disabled={}",
                        cut_group_id.0,
                        cut_id.0,
                        subspace
                            .get_lmb(all_lmbs)
                            .loop_edges
                            .iter()
                            .map(|edge| edge.0)
                            .collect_vec(),
                        variant.disable,
                    );
                    explicit_group_members
                        .entry(group_id)
                        .or_default()
                        .push((signature, description));
                }
                if threshold_candidate
                    .selected_orientation_esurface_id
                    .is_none()
                {
                    debug!(
                        "Leaving graph '{}' cut {:?} threshold {:?} dormant because it has no selected-orientation instance",
                        self.graph.name, cut_edges, threshold_edges,
                    );
                    continue;
                }

                let association = self.threshold_counterterm_association(
                    threshold_candidate,
                    &sandwich,
                    &cut.cut,
                );
                if settings.threshold_subtraction.skip_thresholds_that_are_cuts
                    && self
                        .cuts
                        .iter()
                        .any(|physical_cut| physical_cut.cut == threshold_candidate.cut)
                {
                    continue;
                }
                let configured_variants = Self::configured_threshold_variants(
                    &self.graph.threshold_counterterms,
                    &cut_edges,
                    &threshold_edges,
                );
                let mut variants = Vec::with_capacity(configured_variants.len());
                for (variant, origin) in configured_variants {
                    if !variant.disable
                        && variant
                            .multiplier
                            .as_ref()
                            .is_some_and(|multiplier| multiplier.symmetrize)
                    {
                        unimplemented!(
                            "symmetrized threshold-counterterm multipliers are not implemented (graph '{}', cut {:?}, threshold {:?}, variant '{}')",
                            self.graph.name,
                            cut_edges,
                            threshold_edges,
                            variant.name.as_deref().unwrap_or("default"),
                        );
                    }
                    let name = variant
                        .name
                        .clone()
                        .unwrap_or_else(|| "default".to_string());
                    let context = format!(
                        "graph '{}' cut {:?} threshold {:?} variant '{}'",
                        self.graph.name,
                        cut_edges.iter().map(|edge| edge.0).collect_vec(),
                        threshold_edges.iter().map(|edge| edge.0).collect_vec(),
                        name,
                    );
                    if origin == ThresholdCountertermOrigin::Explicit
                        && variant.parent_lmb.is_none()
                    {
                        return Err(eyre!("{context} requires a full parent_lmb"));
                    }
                    let subspace = self.resolve_threshold_variant_subspace(
                        &variant,
                        legacy_subspace,
                        containing_subgraph,
                        &context,
                    )?;
                    let status = self.classify_threshold_counterterm_association(
                        &association,
                        cut_id,
                        cut,
                        &subspace,
                        all_lmbs,
                        model,
                        runtime_settings,
                        settings,
                    );
                    variants.push(ResolvedThresholdCountertermDraft {
                        name,
                        group_id: variant.group_id,
                        origin,
                        disable: variant.disable,
                        requested_subspace: variant.subspace,
                        requested_parent_lmb: variant.parent_lmb,
                        multiplier: variant.multiplier,
                        subspace,
                        status,
                    });
                }

                let is_eligible = |status: ThresholdCountertermStatus| {
                    status.is_eligible_for_generation(
                        settings.threshold_subtraction.check_esurface_at_generation,
                    )
                };
                let legacy_status = self.classify_threshold_counterterm_association(
                    &association,
                    cut_id,
                    cut,
                    legacy_subspace,
                    all_lmbs,
                    model,
                    runtime_settings,
                    settings,
                );
                if is_eligible(legacy_status)
                    && (variants.len() != 1
                        || variants[0].disable
                        || variants[0].group_id.is_some()
                        || variants[0].multiplier.is_some()
                        || !Self::threshold_subspaces_are_equivalent(
                            &variants[0].subspace,
                            legacy_subspace,
                            all_lmbs,
                        ))
                {
                    legacy_equivalent = false;
                }

                let has_active_variant = variants
                    .iter()
                    .any(|variant| !variant.disable && is_eligible(variant.status));
                if has_active_variant {
                    match side {
                        ThresholdCountertermSide::Amplitude => {
                            unreachable!("amplitude side is not valid in cross-section resolution")
                        }
                        ThresholdCountertermSide::Left => cut_threshold_associations[cut_id]
                            .left
                            .push(association.clone()),
                        ThresholdCountertermSide::Right => cut_threshold_associations[cut_id]
                            .right
                            .push(association.clone()),
                    }
                }

                let raised_threshold_id = raised_threshold_ids[association.esurface_id]
                    .ok_or_else(|| {
                        eyre!(
                            "Graph '{}' threshold E-surface {} is missing from raised-threshold data",
                            self.graph.name,
                            association.esurface_id.0,
                        )
                    })?;
                physical_drafts.push(PhysicalThresholdCountertermDraft {
                    cut_id,
                    cut_group_id,
                    cut_edges: cut_edges.clone(),
                    side,
                    esurface_id: association.esurface_id,
                    threshold_edges,
                    raised_threshold_id,
                    raised_esurface_group: threshold_raised_data.raised_groups[raised_threshold_id]
                        .clone(),
                    variants,
                });
            }
        }

        for configured_cut in &self.graph.threshold_counterterms.cuts {
            if !matched_explicit_cuts.contains(&configured_cut.edges) {
                return Err(eyre!(
                    "Graph '{}' threshold_counterterms cut {:?} does not match any current physical cut",
                    self.graph.name,
                    configured_cut.edges,
                ));
            }
            for configured_threshold in &configured_cut.thresholds {
                if !matched_explicit_thresholds.contains(&(
                    configured_cut.edges.clone(),
                    configured_threshold.edges.clone(),
                )) {
                    return Err(eyre!(
                        "Graph '{}' threshold_counterterms cut {:?} threshold {:?} does not match a topology-discovered association",
                        self.graph.name,
                        configured_cut.edges,
                        configured_threshold.edges,
                    ));
                }
            }
        }

        // An explicit group id is a user assertion about one solve space. Validate every
        // resolved member in its native parent so incompatible cuts can never silently share
        // an SOCP center, including members currently disabled or ineligible for generation.
        let conflicts = explicit_group_members
            .into_iter()
            .filter_map(|(group_id, members)| {
                if members
                    .iter()
                    .all(|(signature, _)| signature == &members[0].0)
                {
                    return None;
                }
                Some(format!(
                    "group_id={group_id} contains incompatible solve subspaces:\n{}",
                    members
                        .into_iter()
                        .map(|(_, description)| description)
                        .join("\n"),
                ))
            })
            .collect_vec();
        if !conflicts.is_empty() {
            return Err(eyre!(
                "Graph '{}' has incompatible explicit threshold groups. A group_id may only refine one automatic signed-cycle subspace. All members of each conflicting group:\n{}",
                self.graph.name,
                conflicts.join("\n"),
            ));
        }

        let mut grouped = BTreeMap::<
            (CutGroupId, RaisedEsurfaceId),
            Vec<PhysicalThresholdCountertermDraft>,
        >::new();
        for draft in physical_drafts {
            grouped
                .entry((draft.cut_group_id, draft.raised_threshold_id))
                .or_default()
                .push(draft);
        }

        let mut aggregated = Vec::new();
        for ((cut_group_id, raised_threshold_id), mut drafts) in grouped {
            let raised_group_is_legacy_equivalent = drafts.iter().all(|draft| {
                let legacy_subspace = match draft.side {
                    ThresholdCountertermSide::Amplitude => return false,
                    ThresholdCountertermSide::Left => {
                        &self.derived_data.subspace_data[cut_group_id].0
                    }
                    ThresholdCountertermSide::Right => {
                        &self.derived_data.subspace_data[cut_group_id].1
                    }
                };
                draft.variants.len() == 1
                    && !draft.variants[0].disable
                    && draft.variants[0].group_id.is_none()
                    && draft.variants[0].multiplier.is_none()
                    && Self::threshold_subspaces_are_equivalent(
                        &draft.variants[0].subspace,
                        legacy_subspace,
                        all_lmbs,
                    )
            });
            let has_left = drafts
                .iter()
                .any(|draft| draft.side == ThresholdCountertermSide::Left);
            let has_right = drafts
                .iter()
                .any(|draft| draft.side == ThresholdCountertermSide::Right);
            if has_left && has_right {
                if raised_group_is_legacy_equivalent {
                    // Preserve the legacy left-side precedence for a raised threshold seen on
                    // both sides of physical members of one cut group.
                    drafts.retain(|draft| draft.side == ThresholdCountertermSide::Left);
                } else {
                    return Err(eyre!(
                        "Graph '{}' cut group {} raised threshold {:?} is inferred on both sides; non-legacy directives must resolve compatibly across raised physical cuts",
                        self.graph.name,
                        cut_group_id.0,
                        threshold_raised_data.raised_groups[raised_threshold_id].esurface_ids,
                    ));
                }
            }

            let representative = drafts.first().expect("grouped drafts are non-empty");
            if drafts.iter().skip(1).any(|draft| {
                draft.variants.len() != representative.variants.len()
                    || draft
                        .variants
                        .iter()
                        .zip(&representative.variants)
                        .any(|(variant, expected)| !variant.is_compatible_with(expected, all_lmbs))
            }) {
                let descriptions = drafts.iter().flat_map(|draft| {
                    draft.variants.iter().map(|variant| format!(
                        "  cut={}, cut_edges={:?}, threshold_edges={:?}, variant='{}', group_id={:?}, parent={:?}, signed_cycles={:?}, disabled={}, multiplier={:?}",
                        draft.cut_id.0,
                        draft.cut_edges,
                        draft.threshold_edges,
                        variant.name,
                        variant.group_id,
                        variant.subspace.get_lmb(all_lmbs).loop_edges,
                        variant.subspace.solve_signature(all_lmbs),
                        variant.disable,
                        variant.multiplier,
                    ))
                }).join("\n");
                return Err(eyre!(
                    "Graph '{}' cut group {} raised threshold {:?} resolves incompatibly across physical cuts. Corresponding counterterm variants of one merged raised residue must have identical directives, group_id and signed solve cycles. All physical members:\n{}",
                    self.graph.name,
                    cut_group_id.0,
                    representative.raised_esurface_group.esurface_ids,
                    descriptions,
                ));
            }
            aggregated.push(drafts);
        }

        aggregated.sort_by_key(|drafts| {
            let representative = &drafts[0];
            (
                representative.cut_group_id.0,
                match representative.side {
                    ThresholdCountertermSide::Amplitude => 2,
                    ThresholdCountertermSide::Left => 0,
                    ThresholdCountertermSide::Right => 1,
                },
                representative
                    .raised_esurface_group
                    .esurface_ids
                    .iter()
                    .map(|id| id.0)
                    .collect_vec(),
            )
        });

        let mut variants =
            TiVec::<ThresholdCountertermVariantId, ResolvedThresholdCountertermVariant>::new();
        let mut cut_groups = ti_vec![
            ResolvedCutGroupThresholdCounterterms::default();
            self.derived_data.cut_group_data.cut_groups.len()
        ];
        let is_eligible = |status: ThresholdCountertermStatus| {
            status.is_eligible_for_generation(
                settings.threshold_subtraction.check_esurface_at_generation,
            )
        };

        for drafts in aggregated {
            let representative = &drafts[0];
            let mut active_variant_count = 0;
            for variant_index in 0..representative.variants.len() {
                let representative_variant = &representative.variants[variant_index];
                if representative_variant.disable {
                    continue;
                }
                let eligible_cut_ids = drafts
                    .iter()
                    .filter(|draft| is_eligible(draft.variants[variant_index].status))
                    .map(|draft| draft.cut_id)
                    .sorted_by_key(|cut_id| cut_id.0)
                    .dedup()
                    .collect_vec();
                if eligible_cut_ids.is_empty() {
                    continue;
                }

                active_variant_count += 1;
                let associations = drafts
                    .iter()
                    .map(|draft| {
                        let resolved_variant = &draft.variants[variant_index];
                        let legacy_subspace = match draft.side {
                            ThresholdCountertermSide::Amplitude => unreachable!(
                                "amplitude side is not valid in cross-section resolution"
                            ),
                            ThresholdCountertermSide::Left => {
                                &self.derived_data.subspace_data[draft.cut_group_id].0
                            }
                            ThresholdCountertermSide::Right => {
                                &self.derived_data.subspace_data[draft.cut_group_id].1
                            }
                        };
                        ResolvedThresholdCountertermAssociation {
                            cut_id: Some(draft.cut_id),
                            cut_edges: draft.cut_edges.clone(),
                            threshold_edges: draft.threshold_edges.clone(),
                            esurface_id: draft.esurface_id,
                            subspace: resolved_variant.subspace.clone(),
                            requires_explicit_parent_lmb: resolved_variant
                                .subspace
                                .parent_lmb_index()
                                != legacy_subspace.parent_lmb_index(),
                            eligible: is_eligible(resolved_variant.status),
                            origin: resolved_variant.origin,
                        }
                    })
                    .sorted_by_key(|association| {
                        (
                            association.cut_id.map(|cut_id| cut_id.0),
                            association.esurface_id.0,
                            association.threshold_edges.clone(),
                        )
                    })
                    .collect();
                let parent_lmb_index = representative_variant.subspace.parent_lmb_index();
                let variant_id = ThresholdCountertermVariantId::from(variants.len());
                variants.push(ResolvedThresholdCountertermVariant {
                    name: representative_variant.name.clone(),
                    group_id: representative_variant.group_id,
                    cut_group_id: Some(representative.cut_group_id),
                    associations,
                    side: representative.side,
                    threshold_esurface_ids: representative
                        .raised_esurface_group
                        .esurface_ids
                        .clone(),
                    raised_esurface_group: representative.raised_esurface_group.clone(),
                    requested_subspace: representative_variant.requested_subspace.clone(),
                    requested_parent_lmb: representative_variant.requested_parent_lmb.clone(),
                    resolved_parent_lmb: all_lmbs[parent_lmb_index].loop_edges.clone().into(),
                    subspace: representative_variant.subspace.clone(),
                    subspace_loop_count: representative_variant.subspace.loopcount(),
                    multiplier: representative_variant.multiplier.clone(),
                });
                match representative.side {
                    ThresholdCountertermSide::Amplitude => {
                        unreachable!("amplitude side is not valid in cross-section resolution")
                    }
                    ThresholdCountertermSide::Left => cut_groups[representative.cut_group_id]
                        .left
                        .push(variant_id),
                    ThresholdCountertermSide::Right => cut_groups[representative.cut_group_id]
                        .right
                        .push(variant_id),
                }
            }

            let legacy_subspace = match representative.side {
                ThresholdCountertermSide::Amplitude => {
                    unreachable!("amplitude side is not valid in cross-section resolution")
                }
                ThresholdCountertermSide::Left => {
                    &self.derived_data.subspace_data[representative.cut_group_id].0
                }
                ThresholdCountertermSide::Right => {
                    &self.derived_data.subspace_data[representative.cut_group_id].1
                }
            };
            if active_variant_count > 0
                && (active_variant_count != 1
                    || representative.variants.len() != 1
                    || representative.variants[0].multiplier.is_some()
                    || representative.variants[0].group_id.is_some()
                    || !Self::threshold_subspaces_are_equivalent(
                        &representative.variants[0].subspace,
                        legacy_subspace,
                        all_lmbs,
                    ))
            {
                legacy_equivalent = false;
            }
        }

        for (cut_group_id, cut_group) in cut_groups.iter_enumerated() {
            for (left_variant_id, right_variant_id) in cut_group.iter_pairs() {
                let left = &variants[left_variant_id];
                let right = &variants[right_variant_id];
                if !left.subspace.is_mergable_with(&right.subspace, all_lmbs) {
                    return Err(eyre!(
                        "Graph '{}' cut group {} has incompatible iterated threshold projections. Left variant '{}': associations={:?}, parent={:?}, signed_cycles={:?}; right variant '{}': associations={:?}, parent={:?}, signed_cycles={:?}. Each selected cycle must have zero response on the other variant's defining edges",
                        self.graph.name,
                        cut_group_id.0,
                        left.name,
                        left.associations
                            .iter()
                            .map(|association| (
                                &association.cut_edges,
                                &association.threshold_edges
                            ))
                            .collect_vec(),
                        left.subspace.get_lmb(all_lmbs).loop_edges,
                        left.subspace.solve_signature(all_lmbs),
                        right.name,
                        right
                            .associations
                            .iter()
                            .map(|association| (
                                &association.cut_edges,
                                &association.threshold_edges
                            ))
                            .collect_vec(),
                        right.subspace.get_lmb(all_lmbs).loop_edges,
                        right.subspace.solve_signature(all_lmbs),
                    ));
                }
            }
        }

        Ok((
            ResolvedThresholdCounterterms {
                legacy_equivalent,
                variants,
                cross_section_cut_groups: cut_groups,
            },
            cut_threshold_associations,
        ))
    }

    fn build_threshold_counterterm(
        &mut self,
        model: &Model,
        settings: &GenerationSettings,
        runtime_settings: &RuntimeSettings,
        all_possible_thresholds: &[TopologicalThresholdCandidate],
        vakint: &Vakint,
        cff_energy_degree_bound_reports: &Mutex<Vec<CffEnergyDegreeBoundReport>>,
    ) -> Result<()> {
        let (resolved, cut_threshold_associations) = self
            .resolve_threshold_counterterm_directives(
                model,
                settings,
                runtime_settings,
                all_possible_thresholds,
            )?;

        let mut cut_structure = vec![];
        let expression_esurfaces = &self
            .derived_data
            .global_cff_expression
            .as_ref()
            .expect("global_cff_expression should have been created")
            .expression
            .surfaces
            .esurface_cache;
        for (cut_group_id, cut_group) in self
            .derived_data
            .cut_group_data
            .cut_groups
            .iter_enumerated()
        {
            let resolved_group = &resolved.cross_section_cut_groups[cut_group_id];
            let lu_cut_selection = cut_group.lu_cut_selection(&self.graph, &self.cuts);
            let cutkosky_cut_union = cut_group
                .cuts
                .iter()
                .map(|cut_id| self.cuts[*cut_id].cut.as_subgraph())
                .reduce(|a, b| a.union(&b))
                .unwrap_or(self.graph.empty_subgraph());

            let add_threshold_group_to_union =
                |base: SuBitGraph, raised_group: &RaisedEsurfaceGroup| {
                    let representative_esurface =
                        &expression_esurfaces[raised_group.esurface_ids[0]];

                    representative_esurface
                        .energies
                        .iter()
                        .map(|edge_id| self.graph.get_edge_subgraph(*edge_id))
                        .fold(base, |acc, subgraph| acc.union(&subgraph))
                };

            for &variant_id in &resolved_group.left {
                let raised_esurface_group = &resolved.variants[variant_id].raised_esurface_group;
                let esurface_cut_union =
                    add_threshold_group_to_union(cutkosky_cut_union.clone(), raised_esurface_group);

                cut_structure.push(CutSet {
                    residue_selector: ResidueSelector {
                        lu: Some(lu_cut_selection.clone()),
                        left_th_cut: Some(raised_esurface_group.clone()),
                        right_th_cut: None,
                    },
                    union: esurface_cut_union,
                    canonicalize_external_shifts: false,
                });
            }

            for &variant_id in &resolved_group.right {
                let raised_esurface_group = &resolved.variants[variant_id].raised_esurface_group;
                let esurface_cut_union =
                    add_threshold_group_to_union(cutkosky_cut_union.clone(), raised_esurface_group);

                cut_structure.push(CutSet {
                    residue_selector: ResidueSelector {
                        lu: Some(lu_cut_selection.clone()),
                        left_th_cut: None,
                        right_th_cut: Some(raised_esurface_group.clone()),
                    },
                    union: esurface_cut_union,
                    canonicalize_external_shifts: false,
                });
            }

            for (left_variant_id, right_variant_id) in resolved_group.iter_pairs() {
                let left_raised_esurface_group =
                    &resolved.variants[left_variant_id].raised_esurface_group;
                let right_raised_esurface_group =
                    &resolved.variants[right_variant_id].raised_esurface_group;
                let esurface_cut_union = add_threshold_group_to_union(
                    add_threshold_group_to_union(
                        cutkosky_cut_union.clone(),
                        left_raised_esurface_group,
                    ),
                    right_raised_esurface_group,
                );

                cut_structure.push(CutSet {
                    residue_selector: ResidueSelector {
                        lu: Some(lu_cut_selection.clone()),
                        left_th_cut: Some(left_raised_esurface_group.clone()),
                        right_th_cut: Some(right_raised_esurface_group.clone()),
                    },
                    union: esurface_cut_union,
                    canonicalize_external_shifts: false,
                });
            }
        }

        let cut_structure = CutStructure {
            cuts: cut_structure,
        };

        let cff_options = self.graph.production_cff_3d_expression_options(settings)?;
        let production_expression = self
            .derived_data
            .global_cff_expression
            .as_ref()
            .expect("global_cff_expression should have been created");

        let mut threshold_counterterms = settings
            .uv
            .orchestrator
            .parametric_integrands(
                &mut self.graph,
                cut_structure,
                vakint,
                OrientationProjection::exact_expression(
                    production_expression,
                    &cff_options,
                    &settings.orientation_pattern,
                    settings.explicit_orientation_sum_only,
                )
                .with_vacuum_subtraction(settings.medium.vacuum_subtraction)
                .with_energy_degree_bound_reports(cff_energy_degree_bound_reports),
                settings,
            )?
            .into_iter();

        let mut result = TiVec::<CutGroupId, LUCounterTermData>::new();
        for (cut_group_id, _) in self
            .derived_data
            .cut_group_data
            .cut_groups
            .iter_enumerated()
        {
            let resolved_group = &resolved.cross_section_cut_groups[cut_group_id];
            let lu_prefactor = self.derived_data.cut_group_data.cut_groups[cut_group_id]
                .lu_prefactor(&self.graph, &self.cuts)?;
            let left_variant_ids: TiVec<LeftThresholdId, ThresholdCountertermVariantId> =
                resolved_group.left.iter().copied().collect();
            let right_variant_ids: TiVec<RightThresholdId, ThresholdCountertermVariantId> =
                resolved_group.right.iter().copied().collect();
            let left_thresholds = left_variant_ids
                .iter()
                .map(|variant_id| resolved.variants[*variant_id].raised_esurface_group.clone())
                .collect();
            let right_thresholds = right_variant_ids
                .iter()
                .map(|variant_id| resolved.variants[*variant_id].raised_esurface_group.clone())
                .collect();
            let left_subspaces = left_variant_ids
                .iter()
                .map(|variant_id| resolved.variants[*variant_id].subspace.clone())
                .collect();
            let right_subspaces = right_variant_ids
                .iter()
                .map(|variant_id| resolved.variants[*variant_id].subspace.clone())
                .collect();

            let mut left_atoms = TiVec::<LeftThresholdId, _>::new();
            let mut right_atoms = TiVec::<RightThresholdId, _>::new();
            let mut iterated_atoms = vec![];

            for variant_id in &left_variant_ids {
                let variant = &resolved.variants[*variant_id];
                let rstar_pow = Self::threshold_rstar_power(variant.subspace_loop_count, false);
                let counterterm = threshold_counterterms.next().ok_or_else(|| {
                    eyre!(
                        "Threshold orchestrator returned too few single-left integrands for graph '{}' cut group {}",
                        self.graph.name,
                        cut_group_id.0,
                    )
                })?;
                left_atoms.push(counterterm.map(|x| x * &lu_prefactor * &rstar_pow));
            }

            for variant_id in &right_variant_ids {
                let variant = &resolved.variants[*variant_id];
                let rstar_pow = Self::threshold_rstar_power(variant.subspace_loop_count, true);
                let counterterm = threshold_counterterms.next().ok_or_else(|| {
                    eyre!(
                        "Threshold orchestrator returned too few single-right integrands for graph '{}' cut group {}",
                        self.graph.name,
                        cut_group_id.0,
                    )
                })?;
                right_atoms.push(counterterm.map(|x| x * &lu_prefactor * &rstar_pow));
            }

            for (left_variant_id, right_variant_id) in resolved_group.iter_pairs() {
                let left_variant = &resolved.variants[left_variant_id];
                let right_variant = &resolved.variants[right_variant_id];
                let left_rstar_pow =
                    Self::threshold_rstar_power(left_variant.subspace_loop_count, false);
                let right_rstar_pow =
                    Self::threshold_rstar_power(right_variant.subspace_loop_count, true);
                let counterterm = threshold_counterterms.next().ok_or_else(|| {
                    eyre!(
                        "Threshold orchestrator returned too few iterated integrands for graph '{}' cut group {}",
                        self.graph.name,
                        cut_group_id.0,
                    )
                })?;
                iterated_atoms.push(
                    counterterm.map(|x| x * &lu_prefactor * &left_rstar_pow * &right_rstar_pow),
                );
            }

            let iterated_collection =
                IteratedCtCollection::new(iterated_atoms, left_atoms.len(), right_atoms.len());

            let counterterm_data = LUCounterTermData {
                left_thresholds,
                right_thresholds,
                left_variant_ids,
                right_variant_ids,
                left_subspaces,
                right_subspaces,
                left_atoms,
                right_atoms,
                iterated: iterated_collection,
            };
            result.push(counterterm_data);
        }
        if threshold_counterterms.next().is_some() {
            return Err(eyre!(
                "Threshold orchestrator returned too many integrands for graph '{}'",
                self.graph.name,
            ));
        }

        let max_dual_size =
            max_dual_size_for_cut_cff_indices(result.iter().flat_map(|counterterm_data| {
                counterterm_data
                    .left_atoms
                    .iter()
                    .chain(counterterm_data.right_atoms.iter())
                    .chain(counterterm_data.iterated.iter())
                    .flat_map(|integrands| integrands.integrands.iter().map(|(index, _)| index))
            }));
        self.graph.param_builder.initialize_duals(max_dual_size);

        self.derived_data.cut_threshold_associations = cut_threshold_associations;
        self.derived_data.threshold_counterterms = result;
        self.derived_data.resolved_threshold_counterterms = Some(resolved);
        Ok(())
    }

    fn build_subspace_data(&mut self) -> Result<()> {
        let all_lmbs = self.derived_data.lmbs.as_ref().unwrap();

        let subspace_data = self
            .derived_data
            .cut_group_data
            .cut_groups
            .iter_enumerated()
            .map(|(cut_group_id, _)| {
                let valid_subspace_lmbs =
                    self.cut_group_compatible_parent_lmbs(cut_group_id)?;
                let (smallest_left_subgraph, smallest_right_subgraph) =
                    self.cut_group_side_subgraphs(cut_group_id)?;

                let mut possible_subspaces = Vec::new();
                let mut rejected_lmbs = Vec::new();
                for lmb_index in valid_subspace_lmbs {
                    let left = SubspaceData::new_with_user_selected_lmb(
                        smallest_left_subgraph.clone(),
                        lmb_index,
                        &self.graph,
                        all_lmbs,
                    );
                    let right = SubspaceData::new_with_user_selected_lmb(
                        smallest_right_subgraph.clone(),
                        lmb_index,
                        &self.graph,
                        all_lmbs,
                    );

                    match (left, right) {
                        (Ok(left), Ok(right)) if left.is_mergable_with(&right, all_lmbs) => {
                            possible_subspaces.push((left, right));
                        }
                        (Ok(left), Ok(right)) => rejected_lmbs.push(format!(
                            "LMB {} produced non-disjoint subspaces: left={:?}, right={:?}",
                            usize::from(lmb_index),
                            left.iter_lmb_indices().collect_vec(),
                            right.iter_lmb_indices().collect_vec(),
                        )),
                        (left, right) => rejected_lmbs.push(format!(
                            "LMB {}: left={}, right={}",
                            usize::from(lmb_index),
                            left.err()
                                .map(|error| format!("{error:#}"))
                                .unwrap_or_else(|| "compatible".to_string()),
                            right
                                .err()
                                .map(|error| format!("{error:#}"))
                                .unwrap_or_else(|| "compatible".to_string()),
                        )),
                    }
                }

                possible_subspaces.sort_by_key(|(left, right)| {
                    (
                        left.iter_basis_edges(all_lmbs).collect_vec(),
                        right.iter_basis_edges(all_lmbs).collect_vec(),
                    )
                });

                possible_subspaces.first().cloned().ok_or_else(|| {
                    eyre!(
                        "No topology-compatible parent LMB found for graph '{}' cut group {}. Rejections:\n{}",
                        self.graph.name,
                        cut_group_id.0,
                        rejected_lmbs.join("\n"),
                    )
                })
            })
            .collect::<Result<_>>()?;

        let _: () = self.derived_data.subspace_data = subspace_data;
        Ok(())
    }

    fn generate_term_for_graph(
        &self,
        _model: &Model,
        settings: &GlobalSettings,
    ) -> Result<(CrossSectionGraphTerm, GraphGenerationStats)> {
        CrossSectionGraphTerm::from_cross_section_graph(self, settings)
    }
}

#[derive(Clone, Encode, Decode)]
#[trait_decode(trait = GammaLoopContext)]
pub struct CrossSectionDerivedData {
    /// Observable labels of the physical states represented by covariant
    /// Goldstone/ghost cuts; explicit unphysical diagnostic states stay intact.
    pub covariant_cut_representatives: BTreeMap<isize, isize>,
    pub orientations: Option<TiVec<OrientationID, EdgeVec<Orientation>>>,
    pub cut_paramatric_integrand: TiVec<CutGroupId, ParametricIntegrands>,
    pub global_cff_expression: Option<
        GeneratedThreeDExpression<crate::cff::esurface::Esurface, crate::cff::hsurface::Hsurface>,
    >,
    pub lmbs: Option<TiVec<LmbIndex, LoopMomentumBasis>>,
    pub multi_channeling_setup: Option<LmbMultiChannelingSetup>,
    pub threshold_counterterms: TiVec<CutGroupId, LUCounterTermData>,
    pub resolved_threshold_counterterms: Option<ResolvedThresholdCounterterms>,
    pub topological_threshold_esurfaces: TiVec<TopologicalThresholdId, Esurface>,
    /// Graph-level inventory of topology-discovered threshold E-surfaces that have an instance
    /// in the selected generation orientations. Orientation-excluded topology identities remain
    /// available while resolving directives, but have no CFF E-surface ID to store here.
    pub threshold_candidate_esurface_ids: Vec<EsurfaceID>,

    /// Exact generated left/right threshold associations for each physical cut. Runtime
    /// evaluators aggregate these into cut groups, but display and model reclassification
    /// must retain the physical-cut-relative eligibility information.
    pub cut_threshold_associations: TiVec<CutId, CutThresholdCountertermAssociations>,
    pub subspace_data: TiVec<CutGroupId, (SubspaceData, SubspaceData)>,
    pub cut_group_data: CutGroupData,
}

/// Partition and residue-derivative metadata for all Cutkosky cuts, including order-one cuts.
#[derive(Clone, Encode, Decode, Debug)]
#[trait_decode(trait = GammaLoopContext)]
pub struct CutGroupData {
    pub cut_groups: TiVec<CutGroupId, CutGroup>,
    pub dual_shapes: Vec<Vec<Vec<usize>>>,
    pub pass_two_evaluators: Vec<GenericEvaluator>,
}

/// Physical cuts that share one normalized E-surface residue, whether ordinary or raised.
#[derive(Clone, Encode, Decode, Debug)]
#[trait_decode(trait = GammaLoopContext)]
pub struct CutGroup {
    pub cuts: Vec<CutId>,
    pub related_esurface_group: RaisedEsurfaceGroup,
}

impl CutGroup {
    pub(crate) fn lu_prefactor(
        &self,
        graph: &Graph,
        cuts: &TiVec<CutId, CrossSectionCut>,
    ) -> Result<Atom> {
        let loop_number =
            graph.cyclotomatic_number(&graph.full_filter()) - graph.initial_state_cut.nedges(graph);

        let loop_3 = loop_number as i64 * 3;
        // The inverse UFO vertices already supply the Hermitian-partner
        // tensors/couplings. Marking the RHS gives (-1)^(V_R+I_R), and its
        // anti-causal virtual contours give (-1)^L_R: together (-1)^C_R.
        // Retain cut hairs when counting components, including a contact RHS.
        let component_parities = self
            .cuts
            .iter()
            .map(|cut| {
                graph
                    .underlying
                    .count_connected_components(&cuts[*cut].right)
                    % 2
            })
            .unique()
            .collect_vec();
        let [parity] = component_parities.as_slice() else {
            return Err(eyre!(
                "Graph '{}' has incompatible RHS marking signs in raised cut group {:?}",
                graph.name,
                self.cuts
            ));
        };
        // CFF retains the cut propagator numerators, including their i factors.
        // The residue conversion supplies 2πi, making the complete physical
        // cut i/(q²-m²+i0) -> 2π δ_+(q²-m²) real. For connected RHS, this
        // combines with marking to -2πi, independent of cut multiplicity.
        let energy_conservation_delta_factor =
            Atom::num(2) * Atom::var(GS.pi) * Atom::i() * Atom::num(-1).pow(*parity as i64);

        crate::debug_tags!(#generation, #normalization, #lu, #graph, #summary;
            stage = "cross_section_lu_prefactor",
            graph = %graph.name,
            loop_number,
            loop_3,
            "Cross-section LU prefactor normalization"
        );

        let tstar = Atom::var(GS.rescale_star);
        let tsrat_pow = tstar.pow(loop_3);
        let hfunction = Atom::var(GS.hfunction_lu_cut);
        Ok(tsrat_pow * hfunction * energy_conservation_delta_factor)
    }

    pub(crate) fn lu_cut_selection(
        &self,
        graph: &Graph,
        cuts: &TiVec<CutId, CrossSectionCut>,
    ) -> LuCutSelection {
        let raised_edge_groups = graph.get_raised_edge_groups();
        let cut_edge_alternatives = self
            .cuts
            .iter()
            .map(|cut_id| {
                let cut_edges = graph
                    .iter_edges_of(&cuts[*cut_id].cut.as_subgraph())
                    .map(|(_, edge_id, _)| edge_id)
                    .collect_vec();
                normalize_cut_edge_support_with_raised_edge_groups(&cut_edges, &raised_edge_groups)
            })
            .collect();
        LuCutSelection {
            raised_group: self.related_esurface_group.clone(),
            cut_edge_alternatives,
        }
    }
}

impl Default for CutGroupData {
    fn default() -> Self {
        Self::new()
    }
}

impl CutGroupData {
    pub fn new() -> Self {
        CutGroupData {
            cut_groups: TiVec::new(),
            dual_shapes: vec![],
            pass_two_evaluators: vec![],
        }
    }

    pub fn new_from_esurface(
        raised_esurface_data: &RaisedEsurfaceData,
        cut_esurface_map: &TiVec<CutId, EsurfaceID>,
        evaluator_settings: &EvaluatorSettings,
    ) -> (Self, GraphGenerationStats) {
        let mut stats = GraphGenerationStats::default();
        let mut groups = TiVec::new();

        for (_raised_esurface_id, raised_esurface_group) in
            raised_esurface_data.raised_groups.iter_enumerated()
        {
            let cuts = cut_esurface_map
                .iter_enumerated()
                .filter_map(|(cut_id, esurface_id)| {
                    raised_esurface_group
                        .esurface_ids
                        .contains(esurface_id)
                        .then_some(cut_id)
                })
                .collect::<Vec<_>>();

            if cuts.is_empty() {
                continue;
            }

            groups.push(CutGroup {
                cuts,
                related_esurface_group: raised_esurface_group.clone(),
            });
        }

        let global_max_occurence = groups
            .iter()
            .map(|group| group.related_esurface_group.max_occurence)
            .max()
            .unwrap_or_else(|| {
                println!("corrupted groups");
                panic!();
            });

        let dual_shapes = (1..global_max_occurence)
            .map(simple_n_deriv_shape)
            .collect();

        let pass_two_evaluators = (1..=global_max_occurence)
            .map(|i| {
                let evaluator_started = std::time::Instant::now();
                let evaluator = build_derivative_structure(i as u8, -1, evaluator_settings);
                stats.evaluator_symbolica_time += evaluator_started.elapsed();
                stats.evaluator_count += 1;
                evaluator
            })
            .collect();

        (
            Self {
                cut_groups: groups,
                dual_shapes,
                pass_two_evaluators,
            },
            stats,
        )
    }
}

impl CrossSectionDerivedData {
    fn new_empty() -> Self {
        Self {
            covariant_cut_representatives: BTreeMap::new(),
            orientations: None,
            global_cff_expression: None,
            cut_paramatric_integrand: TiVec::new(),
            lmbs: None,
            multi_channeling_setup: None,
            threshold_counterterms: TiVec::new(),
            resolved_threshold_counterterms: None,
            topological_threshold_esurfaces: TiVec::new(),
            threshold_candidate_esurface_ids: Vec::new(),

            cut_threshold_associations: TiVec::new(),
            subspace_data: TiVec::new(),
            cut_group_data: CutGroupData::new(),
        }
    }
}

pub(crate) fn build_derivative_structure_atom(
    singularity_order: u8,
    laurent_coefficient: i8,
) -> Atom {
    assert!(
        laurent_coefficient <= -1,
        "only laurent coefficients up to -1 are supported"
    );

    assert!(
        singularity_order >= 1,
        "eta order must be at least 1, got {singularity_order}"
    );
    assert!(
        singularity_order >= -laurent_coefficient as u8,
        "eta order must be at least the negative of the laurent coefficient, got {singularity_order} for laurent coefficient {laurent_coefficient}"
    );

    let order = singularity_order as i32;
    let laurent_coefficient = laurent_coefficient as i32;
    let f = symbol!("f");

    let expansion = function!(GS.eta, GS.rescale)
        .series(GS.rescale, Atom::var(GS.rescale_star), (order, 1))
        .unwrap()
        .to_atom()
        .replace(function!(GS.eta, GS.rescale_star))
        .min_level(0)
        .max_level(Some(0))
        .with(0);

    let mut expression_to_derive = function!(f, GS.rescale)
        * expansion.pow(-order)
        * (GS.rescale - GS.rescale_star).pow(order);

    for _ in 1..=(order + laurent_coefficient) {
        expression_to_derive = expression_to_derive.derivative(GS.rescale);
    }

    expression_to_derive = expression_to_derive
        .replace(GS.rescale - GS.rescale_star)
        .with(parse!("delta_t"));

    let polynomial_in_delta_t = expression_to_derive
        .series(symbol!("delta_t"), Atom::num(0), (0, 1))
        .unwrap();

    let factorial_prefactor = (2..=(order + laurent_coefficient)).product::<i32>();
    debug!("factorial prefactor: {}", factorial_prefactor);
    let mut expression_to_derive = polynomial_in_delta_t.to_atom() / Atom::num(factorial_prefactor);

    expression_to_derive = expression_to_derive
        .replace(GS.rescale)
        .with(GS.rescale_star);

    expression_to_derive
}

pub(crate) fn build_derivative_structure(
    singularity_order: u8,
    laurent_coefficient: i8,
    evaluator_settings: &EvaluatorSettings,
) -> GenericEvaluator {
    let expression_to_derive =
        build_derivative_structure_atom(singularity_order, laurent_coefficient);

    let params = params_for_derivative_order(singularity_order);

    GenericEvaluator::new_from_raw_params(
        [expression_to_derive],
        &params,
        &FunctionMap::default(),
        vec![],
        evaluator_settings.optimization_settings(),
        None,
        evaluator_settings,
    )
    .unwrap()
    .into_eager_only()
}

fn ordered_f_derivative_params(
    base: Atom,
    derivative_vars: &[Symbol],
    derivative_shape: Option<&Vec<Vec<usize>>>,
) -> Vec<Atom> {
    derivative_shape
        .cloned()
        .unwrap_or_else(|| vec![vec![0; derivative_vars.len()]])
        .into_iter()
        .map(|orders| {
            debug_assert_eq!(orders.len(), derivative_vars.len());

            let mut param = base.clone();
            for (derivative_var, order) in derivative_vars.iter().zip(orders) {
                for _ in 0..order {
                    param = param.derivative(*derivative_var);
                }
            }

            param
        })
        .collect()
}

pub(crate) fn params_for_derivative_order(singularity_order: u8) -> Vec<Atom> {
    let f = symbol!("f");

    let f_0 = function!(f, GS.rescale_star);
    let eta_1 = function!(GS.eta, GS.rescale_star).derivative(GS.rescale_star);

    let f_derivative_shape = shape_from_cut_cff_index(&CutCFFIndex {
        left_threshold_order: None,
        right_threshold_order: None,
        lu_cut_order: Some(singularity_order as usize),
    });

    let f_parameters =
        ordered_f_derivative_params(f_0, &[GS.rescale_star], f_derivative_shape.as_ref());
    let mut eta_params = vec![eta_1.clone()];

    for _ in 2..=singularity_order {
        let next_eta = eta_params.last().unwrap().derivative(GS.rescale_star);
        eta_params.push(next_eta);
    }

    let mut result = vec![];
    result.extend(f_parameters);
    result.extend(eta_params);
    result
}

pub(crate) fn params_for_iterated_threshold_ct(
    left_singularit_order: u8,
    right_singularity_order: u8,
) -> Vec<Atom> {
    let f = symbol!("f");

    let cut_cff_index = CutCFFIndex {
        left_threshold_order: Some(left_singularit_order as usize),
        right_threshold_order: Some(right_singularity_order as usize),
        lu_cut_order: None,
    };
    let f_derivative_shape = shape_from_cut_cff_index(&cut_cff_index);

    let f_base = function!(f, GS.radius_star_left, GS.radius_star_right);
    let derivative_vars = match (left_singularit_order > 1, right_singularity_order > 1) {
        (true, true) => vec![GS.radius_star_left, GS.radius_star_right],
        (true, false) => vec![GS.radius_star_left],
        (false, true) => vec![GS.radius_star_right],
        (false, false) => vec![],
    };

    let eta_left_d1 = function!(GS.eta_left, GS.radius_star_left).derivative(GS.radius_star_left);
    let eta_right_d1 =
        function!(GS.eta_right, GS.radius_star_right).derivative(GS.radius_star_right);

    let mut eta_left_params = vec![eta_left_d1.clone()];
    let mut eta_right_params = vec![eta_right_d1.clone()];

    for _ in 2..=left_singularit_order {
        let next_eta_left = eta_left_params
            .last()
            .unwrap()
            .derivative(GS.radius_star_left);

        eta_left_params.push(next_eta_left);
    }

    for _ in 2..=right_singularity_order {
        let next_eta_right = eta_right_params
            .last()
            .unwrap()
            .derivative(GS.radius_star_right);

        eta_right_params.push(next_eta_right);
    }

    let f_product =
        ordered_f_derivative_params(f_base, &derivative_vars, f_derivative_shape.as_ref());

    let mut result = vec![];
    result.extend(f_product);
    result.extend(eta_left_params);
    result.extend(eta_right_params);
    result
}

#[cfg(test)]
mod tests {
    use std::{
        fs,
        path::PathBuf,
        time::{SystemTime, UNIX_EPOCH},
    };

    use symbolica::{atom::AtomCore, function, symbol};

    use crate::{
        cff::CutCFFIndex,
        dot,
        graph::parse::from_dot::IntoGraph,
        initialisation::test_initialise,
        model::InputParamCard,
        processes::cross_section::{EsurfaceID, RaisedEsurfaceData, RaisedEsurfaceGroup},
        settings::global::{GenerationSettings, OrientationPattern},
        utils::{GS, load_generic_model},
    };
    use linnet::half_edge::{
        involution::EdgeIndex,
        subgraph::{OrientedCut, SuBitGraph},
    };

    fn threshold_association(
        cut_boundary_size: usize,
        threshold_boundary_size: usize,
    ) -> super::ThresholdCountertermAssociation {
        super::ThresholdCountertermAssociation {
            topological_threshold_id: super::TopologicalThresholdId(0),
            esurface_id: super::EsurfaceID(0),
            cut_boundary_edges: (0..cut_boundary_size).map(EdgeIndex::from).collect(),
            threshold_boundary_edges: (0..threshold_boundary_size)
                .map(|index| EdgeIndex::from(cut_boundary_size + index))
                .collect(),
            invariant_bound_is_applicable: true,
        }
    }

    #[test]
    fn cut_groups_include_all_physical_cuts_in_a_raised_group() {
        test_initialise().unwrap();
        let raised_esurface_data = RaisedEsurfaceData {
            raised_groups: [
                RaisedEsurfaceGroup {
                    esurface_ids: vec![EsurfaceID(1), EsurfaceID(0)],
                    max_occurence: 1,
                },
                RaisedEsurfaceGroup {
                    esurface_ids: vec![EsurfaceID(2)],
                    max_occurence: 1,
                },
            ]
            .into_iter()
            .collect(),
            pass_two_evaluator: None,
        };
        let cut_esurface_map = [EsurfaceID(0), EsurfaceID(0)].into_iter().collect();

        let (cut_group_data, _) = super::CutGroupData::new_from_esurface(
            &raised_esurface_data,
            &cut_esurface_map,
            &super::EvaluatorSettings::default(),
        );

        let mut cut_partitions: Vec<Vec<super::CutId>> = cut_group_data
            .cut_groups
            .iter()
            .filter(|group| !group.cuts.is_empty())
            .map(|group| {
                let mut cuts = group.cuts.clone();
                cuts.sort_unstable_by_key(|cut| cut.0);
                cuts
            })
            .collect();
        cut_partitions.sort_unstable_by(|left, right| {
            left.iter()
                .map(|cut| cut.0)
                .cmp(right.iter().map(|cut| cut.0))
        });
        assert_eq!(cut_partitions, vec![vec![super::CutId(0), super::CutId(1)]]);
    }

    #[test]
    fn scalar_cut_residue_keeps_the_propagator_i_exactly_once() {
        use crate::graph::{FeynmanGraph, Graph};
        use symbolica::atom::Atom;

        test_initialise().unwrap();
        // This one-line cut is a residue-boundary diagnostic, not a 1->1
        // scattering-rate fixture. It deliberately bypasses the physical
        // process selector's requirement of at least two final particles.
        let graph: Graph = dot!(
            digraph scalar_cut_line {
                node [num=1];
                edge [particle="scalar_0", num="1𝑖"];
                ext [style=invis];
                ext -> left [is_cut=0];
                right -> ext [is_cut=0];
                left -> right;
            }, "scalars"
        )
        .unwrap();
        let mut graph = super::CrossSectionGraph::new(graph);
        graph.cuts = graph
            .graph
            .all_st_cuts_for_cs(
                graph.source_nodes.clone(),
                graph.target_nodes.clone(),
                &graph.graph.get_initial_state_tree(),
            )
            .into_iter()
            .map(|(left, cut, right)| super::CrossSectionCut { cut, left, right })
            .collect();
        assert_eq!(graph.cuts.len(), 1);
        let group = super::CutGroup {
            cuts: vec![super::CutId(0)],
            related_esurface_group: RaisedEsurfaceGroup {
                esurface_ids: vec![EsurfaceID(0)],
                max_occurence: 1,
            },
        };
        let conversion = group
            .lu_prefactor(&graph.graph, &graph.cuts)
            .unwrap()
            .replace(Atom::var(GS.rescale_star))
            .with(Atom::one())
            .replace(Atom::var(GS.hfunction_lu_cut))
            .with(Atom::one());
        let two_pi = Atom::num(2) * Atom::var(GS.pi);
        assert_eq!(conversion, -&two_pi * Atom::i());

        // i/(x+i0)-i/(x-i0)=2*pi*delta(x): the normalized pole
        // coefficient -2*pi*i multiplies the retained numerator i once.
        assert_eq!(&conversion * Atom::i(), two_pi);
        for energy in [1, 2, 5] {
            let q0 = Atom::var(symbol!("scalar_cut_residue_q0"));
            let energy = Atom::num(energy);
            let positive_pole_residue = (Atom::i() / ((&q0 - &energy) * (&q0 + &energy))
                * (&q0 - &energy))
                .replace(q0)
                .with(energy.to_pattern());
            // The delta_+ energy measure is positive pi/E. Inserting an
            // extra cut-line i or reversing the residue sign violates this.
            assert_eq!(
                &conversion * positive_pole_residue,
                Atom::var(GS.pi) / energy
            );
        }
    }

    #[test]
    fn singleton_boundary_classification_uses_generic_mass_hierarchy() {
        use super::ThresholdCountertermStatus;
        use crate::utils::F;

        let association = threshold_association(1, 2);
        let tolerance = F(1.0e-12);

        assert_eq!(
            association.classify_from_invariant_bounds(&F(3.0), &F(4.0), Some(F(3.0)), &tolerance,),
            ThresholdCountertermStatus::ProvenNonExisting,
        );
        assert_eq!(
            association.classify_from_invariant_bounds(&F(4.0), &F(4.0), Some(F(4.0)), &tolerance,),
            ThresholdCountertermStatus::AlwaysPinched,
        );
        assert_eq!(
            association.classify_from_invariant_bounds(&F(5.0), &F(4.0), Some(F(5.0)), &tolerance,),
            ThresholdCountertermStatus::PotentiallyExisting,
        );

        let fixed_association = threshold_association(1, 1);
        assert_eq!(
            fixed_association.classify_from_invariant_bounds(
                &F(3.0),
                &F(4.0),
                Some(F(3.0)),
                &tolerance,
            ),
            ThresholdCountertermStatus::ProvenNonExisting,
        );
        assert_eq!(
            fixed_association.classify_from_invariant_bounds(
                &F(5.0),
                &F(4.0),
                Some(F(5.0)),
                &tolerance,
            ),
            ThresholdCountertermStatus::ProvenNonExisting,
        );
    }

    #[test]
    fn current_model_filter_is_controlled_only_by_generation_setting() {
        use super::ThresholdCountertermStatus;

        assert!(ThresholdCountertermStatus::ProvenNonExisting.is_eligible_for_generation(false));
        assert!(!ThresholdCountertermStatus::ProvenNonExisting.is_eligible_for_generation(true));
        assert!(!ThresholdCountertermStatus::AlwaysPinched.is_eligible_for_generation(false));
        assert!(!ThresholdCountertermStatus::NoRadialDependence.is_eligible_for_generation(false));
        assert!(ThresholdCountertermStatus::PotentiallyExisting.is_eligible_for_generation(true));
    }

    #[test]
    fn multiparticle_boundary_classification_is_conservative() {
        use super::ThresholdCountertermStatus;
        use crate::utils::F;

        let association = threshold_association(2, 2);
        let tolerance = F(1.0e-12);

        assert_eq!(
            association.classify_from_invariant_bounds(&F(4.0), &F(4.0), Some(F(8.0)), &tolerance,),
            ThresholdCountertermStatus::CanBecomePinched,
        );
        assert_eq!(
            association.classify_from_invariant_bounds(&F(3.0), &F(5.0), Some(F(4.0)), &tolerance,),
            ThresholdCountertermStatus::ProvenNonExisting,
        );
        assert_eq!(
            association.classify_from_invariant_bounds(&F(3.0), &F(5.0), Some(F(5.0)), &tolerance,),
            ThresholdCountertermStatus::AlwaysPinched,
        );
        assert_eq!(
            association.classify_from_invariant_bounds(&F(3.0), &F(5.0), Some(F(8.0)), &tolerance,),
            ThresholdCountertermStatus::CanBecomePinched,
        );
        assert_eq!(
            association.classify_from_invariant_bounds(&F(5.0), &F(3.0), Some(F(8.0)), &tolerance,),
            ThresholdCountertermStatus::PotentiallyExisting,
        );
        assert_eq!(
            association.classify_from_invariant_bounds(&F(3.0), &F(5.0), None, &tolerance,),
            ThresholdCountertermStatus::PotentiallyExisting,
        );
    }

    #[test]
    fn invariant_bound_rejects_mixed_or_unresolved_inputs_conservatively() {
        use super::ThresholdCountertermStatus;
        use crate::utils::F;

        let mut association = threshold_association(2, 2);
        let tolerance = F(1.0e-12);
        association.invariant_bound_is_applicable = false;
        assert_eq!(
            association.classify_from_invariant_bounds(&F(3.0), &F(5.0), Some(F(4.0)), &tolerance,),
            ThresholdCountertermStatus::PotentiallyExisting,
        );

        association.invariant_bound_is_applicable = true;
        assert_eq!(
            association.classify_from_invariant_bounds(
                &F(f64::NAN),
                &F(5.0),
                Some(F(4.0)),
                &tolerance,
            ),
            ThresholdCountertermStatus::PotentiallyExisting,
        );

        let masses: linnet::half_edge::involution::EdgeVec<Option<F<f64>>> =
            vec![Some(F(3.0)), None].into();
        assert_eq!(
            super::ThresholdCountertermAssociation::mass_sum(
                &[EdgeIndex::from(0), EdgeIndex::from(1)],
                &masses,
            ),
            None,
        );
    }

    #[test]
    fn singleton_threshold_must_lie_strictly_inside_available_interval() {
        use super::ThresholdCountertermStatus;
        use crate::utils::F;

        let association = threshold_association(2, 1);
        let tolerance = F(1.0e-12);

        assert_eq!(
            association.classify_from_invariant_bounds(&F(5.0), &F(4.0), Some(F(9.0)), &tolerance,),
            ThresholdCountertermStatus::ProvenNonExisting,
        );
        assert_eq!(
            association.classify_from_invariant_bounds(&F(3.0), &F(5.0), Some(F(9.0)), &tolerance,),
            ThresholdCountertermStatus::PotentiallyExisting,
        );
        assert_eq!(
            association.classify_from_invariant_bounds(&F(3.0), &F(9.0), Some(F(9.0)), &tolerance,),
            ThresholdCountertermStatus::AlwaysPinched,
        );
    }

    fn fresh_temp_dir(name: &str) -> PathBuf {
        let unique = SystemTime::now()
            .duration_since(UNIX_EPOCH)
            .unwrap()
            .as_nanos();
        let path = std::env::temp_dir().join(format!(
            "gammalooprs-{name}-{}-{unique}",
            std::process::id()
        ));
        fs::create_dir_all(&path).unwrap();
        path
    }

    #[test]
    fn non_master_cross_section_threshold_metadata_is_rejected_at_import() {
        test_initialise().unwrap();
        let model = load_generic_model("sm");
        let mut master: crate::graph::Graph =
            include_str!("../../../../tests/resources/graphs/ir_safe_thresholds/GL638.dot")
                .into_graph(&model)
                .unwrap();
        master.name = "threshold_group_master".to_string();
        master.group_id = Some(super::GroupId(0));
        master.is_group_master = true;
        let mut member = master.clone();
        member.name = "threshold_group_member".to_string();
        member.is_group_master = false;
        member.threshold_counterterms.autogenerated = false;

        let error = match super::CrossSection::from_graph_list(
            "threshold_group".to_string(),
            vec![master, member],
            &model,
        ) {
            Ok(_) => panic!("non-master threshold metadata must not be interpreted graph-locally"),
            Err(error) => error.to_string(),
        };
        assert!(error.contains("threshold_group_member"));
        assert!(error.contains("threshold_group_master"));
        assert!(error.contains("non-master cross-section graphs is unsupported"));
    }

    #[test]
    fn cross_section_storage_path_stays_inside_process_folder() {
        let temp = fresh_temp_dir("cross-section-storage-path");
        let cross_section = super::CrossSection::new("NLO".to_string());

        assert_eq!(cross_section.storage_path(&temp), temp.join("NLO"));
        fs::remove_dir_all(temp).unwrap();
    }

    #[test]
    fn selected_orientation_keeps_only_selected_threshold_instances_active() {
        use crate::graph::FeynmanGraph;

        std::thread::Builder::new()
            .name("selected-orientation-threshold-test".to_string())
            .stack_size(64 * 1024 * 1024)
            .spawn(|| {
                test_initialise().unwrap();
                let model = load_generic_model("sm");
                let graph =
                    include_str!("../../../../tests/resources/graphs/ir_safe_thresholds/GL638.dot")
                        .into_graph(&model)
                        .unwrap();
                let mut cross_section = super::CrossSectionGraph::new(graph);
                let target_cut_edges =
                    vec![EdgeIndex::from(2), EdgeIndex::from(4), EdgeIndex::from(12)];
                let (left, cut, right) = cross_section
                    .graph
                    .all_st_cuts_for_cs(
                        cross_section.source_nodes.clone(),
                        cross_section.target_nodes.clone(),
                        &cross_section.graph.get_initial_state_tree(),
                    )
                    .into_iter()
                    .find(|(_, candidate, _)| {
                        let mut edges = cross_section
                            .graph
                            .iter_edges_of(candidate)
                            .map(|(_, edge_id, _)| edge_id)
                            .collect::<Vec<_>>();
                        edges.sort();
                        edges == target_cut_edges
                    })
                    .expect("GL638 must contain the target physical cut");
                cross_section
                    .cuts
                    .push(super::CrossSectionCut { cut, left, right });
                let settings = GenerationSettings {
                    orientation_pattern: OrientationPattern::from_user_pattern(
                        "(+,+,+,+,-,-,-,+,-,0,+,0,+,-,+)",
                    )
                    .unwrap(),
                    ..Default::default()
                };
                cross_section.generate_esurface_cuts();
                cross_section.generate_cff(&settings).unwrap();

                let topology_candidates = cross_section
                    .topological_threshold_candidates(&settings)
                    .unwrap();
                let selected_orientation_energy_sets = topology_candidates
                    .iter()
                    .filter_map(|candidate| {
                        candidate
                            .selected_orientation_esurface_id
                            .map(|esurface_id| cross_section.threshold_energy_edges(esurface_id))
                    })
                    .collect::<Vec<_>>();
                for expected in [
                    vec![EdgeIndex::from(7), EdgeIndex::from(8)],
                    vec![EdgeIndex::from(8), EdgeIndex::from(12), EdgeIndex::from(14)],
                    vec![EdgeIndex::from(5), EdgeIndex::from(10)],
                    vec![EdgeIndex::from(5), EdgeIndex::from(12), EdgeIndex::from(13)],
                ] {
                    assert!(selected_orientation_energy_sets.contains(&expected));
                }
                let orientation_excluded_edges = vec![
                    EdgeIndex::from(8),
                    EdgeIndex::from(10),
                    EdgeIndex::from(13),
                    EdgeIndex::from(14),
                ];
                assert!(!selected_orientation_energy_sets.contains(&orientation_excluded_edges));
                assert!(topology_candidates.iter().any(|candidate| {
                    candidate.threshold_esurface.energies == orientation_excluded_edges
                        && candidate.selected_orientation_esurface_id.is_none()
                }));

                // Orientation selection may leave topology-wide candidates dormant, but a surface
                // required by a selected-orientation instance must still exist in the CFF cache.
                cross_section.graph.surface_cache.esurface_cache.raw.clear();
                let error = cross_section
                    .topological_threshold_candidates(&settings)
                    .err()
                    .expect("a corrupted selected CFF cache must be rejected");
                assert!(
                    error
                        .to_string()
                        .contains("has a selected-orientation instance but is missing")
                );
            })
            .unwrap()
            .join()
            .unwrap();
    }

    #[test]
    fn orientation_excluded_explicit_threshold_is_dormant_not_unmatched() {
        std::thread::Builder::new()
            .name("dormant-orientation-threshold-test".to_string())
            .stack_size(64 * 1024 * 1024)
            .spawn(|| {
                use crate::graph::{
                    FeynmanGraph,
                    autogen::Autogen,
                    threshold_counterterms::{
                        ThresholdCountertermCut, ThresholdCountertermMultiplier,
                        ThresholdCountertermSpec, ThresholdCountertermThreshold,
                        ThresholdCountertermVariant,
                    },
                };

                test_initialise().unwrap();
                let model = load_generic_model("sm");
                let graph =
                    include_str!("../../../../tests/resources/graphs/ir_safe_thresholds/GL638.dot")
                        .into_graph(&model)
                        .unwrap();
                let mut cross_section = super::CrossSectionGraph::new(graph);
                let target_cut_edges =
                    vec![EdgeIndex::from(2), EdgeIndex::from(4), EdgeIndex::from(12)];
                let (left, cut, right) = cross_section
                    .graph
                    .all_st_cuts_for_cs(
                        cross_section.source_nodes.clone(),
                        cross_section.target_nodes.clone(),
                        &cross_section.graph.get_initial_state_tree(),
                    )
                    .into_iter()
                    .find(|(_, candidate, _)| {
                        let mut edges = cross_section
                            .graph
                            .iter_edges_of(candidate)
                            .map(|(_, edge_id, _)| edge_id)
                            .collect::<Vec<_>>();
                        edges.sort();
                        edges == target_cut_edges
                    })
                    .expect("GL638 must contain the target physical cut");
                cross_section
                    .cuts
                    .push(super::CrossSectionCut { cut, left, right });

                let dormant_threshold_edges = vec![EdgeIndex::from(5), EdgeIndex::from(10)];
                cross_section.graph.threshold_counterterms =
                    Autogen::explicit(ThresholdCountertermSpec {
                        schema_version: 1,
                        function_map: Default::default(),
                        cuts: vec![ThresholdCountertermCut {
                            edges: target_cut_edges,
                            thresholds: vec![ThresholdCountertermThreshold {
                                edges: dormant_threshold_edges.clone(),
                                counterterms: vec![ThresholdCountertermVariant {
                                    group_id: None,
                                    name: Some("dormant".to_string()),
                                    subspace: None,
                                    parent_lmb: Some(vec![
                                        EdgeIndex(3),
                                        EdgeIndex(4),
                                        EdgeIndex(7),
                                        EdgeIndex(10),
                                    ]),
                                    disable: false,
                                    multiplier: Some(ThresholdCountertermMultiplier {
                                        expression: "1".to_string(),
                                        function_map: Default::default(),
                                        symmetrize: true,
                                        opaque_derivatives: true,
                                    }),
                                }],
                            }],
                        }],
                    });

                // This orientation contains the physical cut and both left thresholds, but excludes the
                // right-side (5,10) threshold. Its explicit directive must remain matched and dormant.
                let settings = GenerationSettings {
                    orientation_pattern: OrientationPattern::from_user_pattern(
                        "(+,+,+,+,-,+,-,+,-,0,+,0,+,-,+)",
                    )
                    .unwrap(),
                    ..Default::default()
                };
                cross_section.generate_esurface_cuts();
                cross_section.generate_cff(&settings).unwrap();
                cross_section.build_lmbs().unwrap();
                cross_section.build_subspace_data().unwrap();
                let topology_candidates = cross_section
                    .topological_threshold_candidates(&settings)
                    .unwrap();
                assert!(topology_candidates.iter().any(|candidate| {
                    candidate.threshold_esurface.energies == dormant_threshold_edges
                        && candidate.selected_orientation_esurface_id.is_none()
                }));

                let (resolved, _) = cross_section
                    .resolve_threshold_counterterm_directives(
                        &model,
                        &settings,
                        &crate::settings::RuntimeSettings::default(),
                        &topology_candidates,
                    )
                    .unwrap();
                assert!(resolved.variants.iter().all(|variant| {
                    variant.name != "dormant"
                        && variant.associations.iter().all(|association| {
                            association.threshold_edges != dormant_threshold_edges
                        })
                }));
                assert!(
                    resolved
                        .cross_section_cut_groups
                        .iter()
                        .any(|cut_group| !cut_group.right.is_empty()),
                    "the orientation retains its other right-side threshold"
                );
            })
            .unwrap()
            .join()
            .unwrap();
    }

    #[test]
    fn cross_section_lmb_cut_edge_exclusion_prefers_massive_then_fermion_then_id() {
        test_initialise().unwrap();
        let mut model = load_generic_model("sm");
        let mut parameters = InputParamCard::from_str(
            include_str!("../../../../assets/models/json/sm/restrict_light_quarks_massless.json")
                .to_owned(),
            "json",
        )
        .unwrap();
        model.simplify(&mut parameters).unwrap();
        // Keep the light fermion massless to distinguish mass priority from spin priority.
        let graph = dot!(
            digraph cut_edge_priority {
                edge [num=1]
                node [num=1]
                A -> B [id=0 particle="g"]
                A -> B [id=1 particle="d"]
                A -> B [id=2 particle="t" mass=1]
                A -> B [id=3 particle="a"]
                A -> B [id=4 mass=0]
            },
            &model
        )
        .unwrap();

        assert_eq!(
            super::CrossSectionGraph::excluded_cut_edge_for_lmb_channel_in(
                &graph,
                &[EdgeIndex::from(0), EdgeIndex::from(1), EdgeIndex::from(2)]
            ),
            Some(EdgeIndex::from(2))
        );
        assert_eq!(
            super::CrossSectionGraph::excluded_cut_edge_for_lmb_channel_in(
                &graph,
                &[EdgeIndex::from(0), EdgeIndex::from(1), EdgeIndex::from(3)]
            ),
            Some(EdgeIndex::from(1))
        );
        assert_eq!(
            super::CrossSectionGraph::excluded_cut_edge_for_lmb_channel_in(
                &graph,
                &[EdgeIndex::from(0), EdgeIndex::from(3)]
            ),
            Some(EdgeIndex::from(0))
        );

        // A structural zero and a DOT-overridden mass are resolved information for current-model
        // invariant bounds. Keep genuinely unresolved and complex masses conservative, but do not
        // let a massless edge disable a proof independent of model parameters.
        assert_eq!(
            graph.underlying[EdgeIndex::from(2)]
                .particle
                .mass_atom()
                .to_string(),
            "1"
        );
        let empty: SuBitGraph = graph.underlying.empty_subgraph();
        let cut = super::CrossSectionCut {
            cut: OrientedCut {
                left: empty.clone(),
                right: empty.clone(),
            },
            left: empty.clone(),
            right: empty,
        };
        let association = super::ThresholdCountertermAssociation {
            topological_threshold_id: super::TopologicalThresholdId(0),
            esurface_id: super::EsurfaceID(0),
            cut_boundary_edges: vec![EdgeIndex::from(4)],
            threshold_boundary_edges: vec![EdgeIndex::from(2)],
            invariant_bound_is_applicable: true,
        };
        assert_eq!(
            association.classify_for_model(
                &graph,
                &cut,
                &crate::model::Model::default(),
                &graph.param_builder,
                &crate::settings::RuntimeSettings::default(),
                1.0e-7,
            ),
            super::ThresholdCountertermStatus::ProvenNonExisting,
        );
    }

    #[test]
    fn iterated_counterterms_use_the_right_threshold_count_as_row_stride() {
        let two_by_three = super::IteratedCtCollection::new((0..6).collect(), 2, 3);
        assert_eq!(
            two_by_three[(
                super::LeftThresholdId::from(0),
                super::RightThresholdId::from(0)
            )],
            0
        );
        assert_eq!(
            two_by_three[(
                super::LeftThresholdId::from(0),
                super::RightThresholdId::from(2)
            )],
            2
        );
        assert_eq!(
            two_by_three[(
                super::LeftThresholdId::from(1),
                super::RightThresholdId::from(0)
            )],
            3
        );
        assert_eq!(
            two_by_three[(
                super::LeftThresholdId::from(1),
                super::RightThresholdId::from(2)
            )],
            5
        );

        let three_by_two = super::IteratedCtCollection::new((0..6).collect(), 3, 2);
        assert_eq!(
            three_by_two[(
                super::LeftThresholdId::from(0),
                super::RightThresholdId::from(1)
            )],
            1
        );
        assert_eq!(
            three_by_two[(
                super::LeftThresholdId::from(1),
                super::RightThresholdId::from(0)
            )],
            2
        );
        assert_eq!(
            three_by_two[(
                super::LeftThresholdId::from(2),
                super::RightThresholdId::from(1)
            )],
            5
        );
    }

    #[test]
    fn directive_resolution_defaults_and_explicit_lists_are_authoritative() {
        use crate::{
            graph::threshold_counterterms::{
                ThresholdCountertermCut, ThresholdCountertermSpec, ThresholdCountertermThreshold,
                ThresholdCountertermVariant,
            },
            processes::threshold_counterterms::ThresholdCountertermOrigin,
        };

        let cut_edges = vec![EdgeIndex::from(2), EdgeIndex::from(4)];
        let threshold_edges = vec![EdgeIndex::from(7), EdgeIndex::from(8)];
        let mut spec = ThresholdCountertermSpec {
            schema_version: 1,
            function_map: Default::default(),
            cuts: vec![ThresholdCountertermCut {
                edges: cut_edges.clone(),
                thresholds: vec![ThresholdCountertermThreshold {
                    edges: threshold_edges.clone(),
                    counterterms: Vec::new(),
                }],
            }],
        };

        let implicit = super::CrossSectionGraph::configured_threshold_variants(
            &spec,
            &cut_edges,
            &threshold_edges,
        );
        assert_eq!(implicit.len(), 1);
        assert_eq!(implicit[0].0.name.as_deref(), Some("default"));
        assert_eq!(implicit[0].1, ThresholdCountertermOrigin::Autogenerated);

        spec.cuts[0].thresholds[0].counterterms = vec![
            ThresholdCountertermVariant {
                group_id: None,
                name: Some("disabled".to_string()),
                subspace: None,
                parent_lmb: Some(vec![EdgeIndex(7), EdgeIndex(8)]),
                disable: true,
                multiplier: None,
            },
            ThresholdCountertermVariant {
                group_id: None,
                name: Some("duplicate".to_string()),
                subspace: Some(vec![EdgeIndex::from(7)]),
                parent_lmb: Some(vec![EdgeIndex(7), EdgeIndex(8)]),
                disable: false,
                multiplier: None,
            },
        ];
        let explicit = super::CrossSectionGraph::configured_threshold_variants(
            &spec,
            &cut_edges,
            &threshold_edges,
        );
        assert_eq!(explicit.len(), 2);
        assert!(explicit[0].0.disable);
        assert_eq!(explicit[1].0.name.as_deref(), Some("duplicate"));
        assert!(
            explicit
                .iter()
                .all(|(_, origin)| *origin == ThresholdCountertermOrigin::Explicit)
        );
    }

    #[test]
    fn explicit_group_conflicts_report_all_members_across_cut_sides() {
        use itertools::Itertools;

        std::thread::Builder::new()
            .name("threshold-group-conflicts".to_string())
            .stack_size(64 * 1024 * 1024)
            .spawn(|| {
                use crate::graph::{
                    FeynmanGraph,
                    autogen::Autogen,
                    threshold_counterterms::{
                        ThresholdCountertermCut, ThresholdCountertermSpec,
                        ThresholdCountertermThreshold, ThresholdCountertermVariant,
                    },
                };

                test_initialise().unwrap();
                let model = load_generic_model("sm");
                let graph =
                    include_str!("../../../../tests/resources/graphs/ir_safe_thresholds/GL297.dot")
                        .into_graph(&model)
                        .unwrap();
                let mut graph = super::CrossSectionGraph::new(graph);
                let declarations = [
                    ([2, 4, 9], [5, 7], [4, 5, 9, 11], 5, "left"),
                    ([2, 6, 7], [3, 9], [3, 6, 7, 11], 3, "right"),
                ];
                graph.cuts = graph
                    .graph
                    .all_st_cuts_for_cs(
                        graph.source_nodes.clone(),
                        graph.target_nodes.clone(),
                        &graph.graph.get_initial_state_tree(),
                    )
                    .into_iter()
                    .filter(|(_, cut, _)| {
                        let edges = graph
                            .graph
                            .iter_edges_of(cut)
                            .map(|(_, edge, _)| edge.0)
                            .sorted()
                            .collect_vec();
                        declarations.iter().any(|(cut, ..)| edges.as_slice() == cut)
                    })
                    .map(|(left, cut, right)| super::CrossSectionCut { cut, left, right })
                    .collect();
                assert_eq!(graph.cuts.len(), 2);
                graph.graph.threshold_counterterms = Autogen::explicit(ThresholdCountertermSpec {
                    schema_version: 1,
                    function_map: Default::default(),
                    cuts: declarations
                        .into_iter()
                        .map(
                            |(cut, threshold, parent, active, side)| ThresholdCountertermCut {
                                edges: cut.into_iter().map(EdgeIndex).collect(),
                                thresholds: vec![ThresholdCountertermThreshold {
                                    edges: threshold.into_iter().map(EdgeIndex).collect(),
                                    counterterms: (0..2)
                                        .map(|group_id| ThresholdCountertermVariant {
                                            name: Some(format!("{side}_group_{group_id}")),
                                            group_id: Some(group_id),
                                            subspace: Some(vec![EdgeIndex(active)]),
                                            parent_lmb: Some(
                                                parent.into_iter().map(EdgeIndex).collect(),
                                            ),
                                            disable: true,
                                            multiplier: None,
                                        })
                                        .collect(),
                                }],
                            },
                        )
                        .collect(),
                });
                let settings = GenerationSettings {
                    orientation_pattern: OrientationPattern::from_user_pattern(
                        "(+,+,-,+,+,+,+,-,0,-,0,-,-,+,+)",
                    )
                    .unwrap(),
                    ..Default::default()
                };
                graph.generate_esurface_cuts();
                graph.generate_cff(&settings).unwrap();
                graph.build_lmbs().unwrap();
                graph.build_subspace_data().unwrap();
                let topology = graph.topological_threshold_candidates(&settings).unwrap();
                let error = graph
                    .resolve_threshold_counterterm_directives(
                        &model,
                        &settings,
                        &crate::settings::RuntimeSettings::default(),
                        &topology,
                    )
                    .expect_err("explicit groups cannot join different cut-side cycles");
                let diagnostic = error.to_string();
                for expected in [
                    "group_id=0",
                    "group_id=1",
                    "left_group_0",
                    "right_group_0",
                    "left_group_1",
                    "right_group_1",
                    "signed_cycles=",
                    "parent=",
                    "cut [2, 4, 9]",
                    "cut [2, 6, 7]",
                    "disabled=true",
                ] {
                    assert!(
                        diagnostic.contains(expected),
                        "missing {expected}: {diagnostic}"
                    );
                }
            })
            .unwrap()
            .join()
            .unwrap();
    }

    #[test]
    fn raised_member_compatibility_uses_directives_and_embedding_not_provenance() {
        use crate::{
            graph::{Graph, LMBext},
            momentum::sample::SubspaceData,
            processes::threshold_counterterms::ThresholdCountertermOrigin,
        };

        test_initialise().unwrap();
        let graph: Graph = dot!(digraph raised_compatibility {
            ext [style=invis]
            edge [num=1 mass=0]
            node [num=1]
            ext->v1:0 [id=0]
            v1->v2 [id=1]
            v2->v1 [id=2]
            v2->v3 [id=3]
            v3->v2 [id=4]
            ext->v3:1 [id=5]
        })
        .unwrap();
        let all_lmbs = graph.generate_loop_momentum_bases();
        let subspaces = all_lmbs
            .iter_enumerated()
            .filter_map(|(lmb_index, _)| {
                SubspaceData::new_from_parent_basis_edges(
                    &[EdgeIndex(1)],
                    &graph.full_filter(),
                    lmb_index,
                    &graph,
                    &all_lmbs,
                )
                .ok()
            })
            .collect::<Vec<_>>();
        let (subspace, raised_subspace) = subspaces
            .iter()
            .enumerate()
            .flat_map(|(left_index, left)| {
                subspaces
                    .iter()
                    .skip(left_index + 1)
                    .map(move |right| (left, right))
            })
            .find(|(left, right)| {
                left.parent_lmb_index() != right.parent_lmb_index()
                    && left.has_equivalent_embedding(right, &all_lmbs)
            })
            .expect("raised fixture must have equivalent distinct parent embeddings");
        let requested = subspace.iter_basis_edges(&all_lmbs).collect::<Vec<_>>();
        let first = super::ResolvedThresholdCountertermDraft {
            name: "v".to_string(),
            group_id: None,
            origin: ThresholdCountertermOrigin::Explicit,
            disable: false,
            requested_subspace: Some(requested),
            requested_parent_lmb: Some(subspace.get_lmb(&all_lmbs).loop_edges.raw.clone()),
            multiplier: None,
            subspace: subspace.clone(),
            status: super::ThresholdCountertermStatus::PotentiallyExisting,
        };
        let mut raised_member = first.clone();
        raised_member.origin = ThresholdCountertermOrigin::Autogenerated;
        raised_member.status = super::ThresholdCountertermStatus::ProvenNonExisting;
        raised_member.requested_subspace = None;
        raised_member.requested_parent_lmb = Some(
            raised_subspace
                .get_lmb(&all_lmbs)
                .loop_edges
                .iter()
                .copied()
                .collect(),
        );
        raised_member.subspace = raised_subspace.clone();

        // Provenance and syntactic requests live on each physical declaration. Raised-member
        // compatibility is semantic, so equivalent resolved embeddings remain compatible.
        assert!(
            super::CrossSectionGraph::threshold_subspaces_are_equivalent(
                subspace,
                raised_subspace,
                &all_lmbs,
            )
        );
        assert!(first.is_compatible_with(&raised_member, &all_lmbs));
        raised_member.name = "different".to_string();
        assert!(!first.is_compatible_with(&raised_member, &all_lmbs));
    }

    #[test]
    fn threshold_radial_power_uses_each_variant_loop_count() {
        let one_loop = super::CrossSectionGraph::threshold_rstar_power(1, false);
        let two_loop = super::CrossSectionGraph::threshold_rstar_power(2, true);

        assert!(
            (one_loop - super::Atom::var(GS.radius_star_left).pow(2))
                .expand()
                .is_zero()
        );
        assert!(
            (two_loop - super::Atom::var(GS.radius_star_right).pow(5))
                .expand()
                .is_zero()
        );
    }

    #[test]
    fn max_dual_size_for_cut_cff_indices_tracks_mixed_threshold_shapes() {
        let cut_cff_indices = [
            CutCFFIndex {
                left_threshold_order: None,
                right_threshold_order: None,
                lu_cut_order: Some(2),
            },
            CutCFFIndex {
                left_threshold_order: Some(2),
                right_threshold_order: Some(2),
                lu_cut_order: Some(2),
            },
        ];

        assert_eq!(
            super::max_dual_size_for_cut_cff_indices(cut_cff_indices.iter()),
            8
        );
    }

    #[test]
    fn single_threshold_params_follow_effective_f_then_eta_contract() {
        let params = super::params_for_derivative_order(3);
        let f = symbol!("f");
        let f_base = function!(f, GS.rescale_star);
        let eta_base = function!(GS.eta, GS.rescale_star);

        let expected = [
            f_base.clone(),
            f_base.clone().derivative(GS.rescale_star),
            f_base
                .derivative(GS.rescale_star)
                .derivative(GS.rescale_star),
            eta_base.clone().derivative(GS.rescale_star),
            eta_base
                .clone()
                .derivative(GS.rescale_star)
                .derivative(GS.rescale_star),
            eta_base
                .derivative(GS.rescale_star)
                .derivative(GS.rescale_star)
                .derivative(GS.rescale_star),
        ]
        .into_iter()
        .map(|atom| atom.to_string())
        .collect::<Vec<_>>();

        let actual = params
            .iter()
            .map(|atom| atom.to_string())
            .collect::<Vec<_>>();

        assert_eq!(actual, expected);
    }

    #[test]
    fn iterated_threshold_f_params_follow_mixed_left_right_shape_order() {
        let params = super::params_for_iterated_threshold_ct(2, 2);
        let f = symbol!("f");
        let base = function!(f, GS.radius_star_left, GS.radius_star_right);
        let expected = [
            base.clone(),
            base.clone().derivative(GS.radius_star_left),
            base.clone().derivative(GS.radius_star_right),
            base.derivative(GS.radius_star_left)
                .derivative(GS.radius_star_right),
        ]
        .into_iter()
        .map(|atom| atom.to_string())
        .collect::<Vec<_>>();
        let actual = params[..expected.len()]
            .iter()
            .map(|atom| atom.to_string())
            .collect::<Vec<_>>();

        assert_eq!(actual, expected);
    }

    #[test]
    fn iterated_threshold_params_append_left_then_right_eta_families() {
        let params = super::params_for_iterated_threshold_ct(2, 2);

        let expected = [
            function!(GS.eta_left, GS.radius_star_left).derivative(GS.radius_star_left),
            function!(GS.eta_left, GS.radius_star_left)
                .derivative(GS.radius_star_left)
                .derivative(GS.radius_star_left),
            function!(GS.eta_right, GS.radius_star_right).derivative(GS.radius_star_right),
            function!(GS.eta_right, GS.radius_star_right)
                .derivative(GS.radius_star_right)
                .derivative(GS.radius_star_right),
        ]
        .into_iter()
        .map(|atom| atom.to_string())
        .collect::<Vec<_>>();

        let actual = params[4..]
            .iter()
            .map(|atom| atom.to_string())
            .collect::<Vec<_>>();

        assert_eq!(actual, expected);
    }

    #[test]
    fn iterated_threshold_helpers_match_signed_nonlinear_pole_distributions() {
        use symbolica::atom::Atom;

        test_initialise().unwrap();
        // For x=r-r*, take f=2+7x+11x²/2 and eta=2x+3x²/2+5x³/6.
        // Direct meromorphic residues of f/eta^n for n=1,2,3 are 1,1,-3/4.
        // The remaining Laurent coefficients are (1/2) for n=2 and
        // (5/16,1/4) for n=3. At r=2,r*=1 and dampers (3,5), their
        // principal-value helpers are 4/3,17/9,-89/216 before r^(1-3L).
        // These fixed scalar jets are not a graph or loop-momentum numerator.
        let residues = [Atom::one(), Atom::one(), Atom::num((-3, 4))];
        let principal_values = [Atom::num((4, 3)), Atom::num((17, 9)), Atom::num((-89, 216))];
        let f_jets = [2_i64, 7, 11];
        let eta_jets = [2_i64, 3, 5];
        let f = symbol!("f");
        for left_order in 1..=3_u8 {
            for right_order in 1..=3_u8 {
                let (_, helper) = super::CrossSectionGraph::iterated_th_prefactor_helper_atoms(
                    left_order,
                    right_order,
                    1,
                    2,
                    true,
                );
                let mut sampled = helper;
                // A separable bivariate test function fixes every mixed jet
                // independently of the production parameter-layout helpers.
                for (left_degree, left_jet) in f_jets.iter().enumerate() {
                    for (right_degree, right_jet) in f_jets.iter().enumerate() {
                        let mut derivative =
                            function!(f, GS.radius_star_left, GS.radius_star_right);
                        for _ in 0..left_degree {
                            derivative = derivative.derivative(GS.radius_star_left);
                        }
                        for _ in 0..right_degree {
                            derivative = derivative.derivative(GS.radius_star_right);
                        }
                        sampled = sampled
                            .replace(derivative)
                            .with(Atom::num(left_jet * right_jet));
                    }
                }
                for (eta, radius_star) in [
                    (GS.eta_left, GS.radius_star_left),
                    (GS.eta_right, GS.radius_star_right),
                ] {
                    let mut derivative = function!(eta, radius_star);
                    for jet in eta_jets {
                        derivative = derivative.derivative(radius_star);
                        sampled = sampled.replace(derivative.clone()).with(Atom::num(jet));
                    }
                }
                for (parameter, value) in [
                    (GS.radius_left, 2),
                    (GS.radius_right, 2),
                    (GS.radius_star_left, 1),
                    (GS.radius_star_right, 1),
                    (GS.uv_damp_plus_left, 3),
                    (GS.uv_damp_plus_right, 3),
                    (GS.uv_damp_minus_left, 5),
                    (GS.uv_damp_minus_right, 5),
                    (GS.pi, 3),
                ] {
                    sampled = sampled.replace(parameter).with(Atom::num(value));
                }
                for (left_integrated, right_integrated) in [(0, 0), (1, 0), (0, 1), (1, 1)] {
                    let actual = sampled
                        .replace(GS.hfunction_left_th)
                        .with(Atom::num(left_integrated))
                        .replace(GS.hfunction_right_th)
                        .with(Atom::num(right_integrated));
                    // Sokhotski-Plemelj fixes the helper signs because runtime
                    // subtracts each single CT and adds their intersection.
                    // The mixed terms have opposite imaginary signs, while
                    // the two delta terms multiply to positive pi².
                    let left = (&principal_values[left_order as usize - 1]
                        - Atom::i()
                            * Atom::num(3 * left_integrated)
                            * &residues[left_order as usize - 1])
                        * Atom::num((1, 4));
                    let right = (&principal_values[right_order as usize - 1]
                        + Atom::i()
                            * Atom::num(3 * right_integrated)
                            * &residues[right_order as usize - 1])
                        * Atom::num((1, 32));
                    assert_eq!(
                        actual,
                        left * right,
                        "orders {left_order}x{right_order}, integrated {left_integrated}/{right_integrated}"
                    );
                }
            }
        }
    }

    #[test]
    fn iterated_threshold_helper_atom_matches_its_left_right_parameter_families() {
        std::thread::Builder::new()
            .name("iterated-threshold-helper-test".to_string())
            .stack_size(16 * 1024 * 1024)
            .spawn(|| {
                test_initialise().unwrap();

                let mut fn_map = super::FunctionMap::new();
                fn_map
                    .add_aliases([(
                        GS.pi.into(),
                        super::Atom::num(
                            super::Rational::try_from(std::f64::consts::PI).unwrap(),
                        ),
                    )])
                    .unwrap();

                for left_order in 1..=3 {
                    for right_order in 1..=3 {
                        let (pieces, legacy) =
                            super::CrossSectionGraph::iterated_th_prefactor_helper_atoms(
                            left_order,
                            right_order,
                            1,
                            1,
                            true,
                        );
                        let params = super::CrossSectionGraph::iterated_th_prefactor_helper_params(
                            left_order,
                            right_order,
                        );

                        super::GenericEvaluator::new_from_raw_params(
                            [
                                legacy,
                                pieces.local_local,
                                pieces.local_integrated,
                                pieces.integrated_local,
                                pieces.integrated_integrated,
                            ],
                            &params,
                            &fn_map,
                            vec![],
                            super::OptimizationSettings::default(),
                            None,
                            &super::EvaluatorSettings::default(),
                        )
                        .unwrap_or_else(|error| {
                            panic!(
                                "iterated {left_order}x{right_order} threshold helper atom and parameter families must agree: {error}"
                            )
                        });
                    }
                }
            })
            .expect("iterated threshold-helper test thread must start")
            .join()
            .expect("iterated threshold-helper test thread must finish successfully");
    }

    #[test]
    fn single_threshold_helper_keeps_raised_poles_in_the_local_piece() {
        let (with_integrated, legacy) =
            super::CrossSectionGraph::single_th_prefactor_helper_atoms(3, 1, false, true);
        let (without_integrated, _) =
            super::CrossSectionGraph::single_th_prefactor_helper_atoms(3, 1, false, false);

        assert!(
            (with_integrated.local.clone() - without_integrated.local)
                .expand()
                .is_zero()
        );
        assert!(without_integrated.integrated.is_zero());
        assert!(!with_integrated.integrated.is_zero());
        assert!(
            (legacy - with_integrated.local - with_integrated.integrated)
                .expand()
                .is_zero()
        );
    }

    #[test]
    fn iterated_threshold_helper_pieces_sum_to_the_legacy_product() {
        let (left, _) =
            super::CrossSectionGraph::single_th_prefactor_helper_atoms(3, 1, false, true);
        let left_local = left
            .local
            .replace(GS.eta)
            .with(super::Atom::var(GS.eta_left));
        let left_integrated = left
            .integrated
            .replace(GS.eta)
            .with(super::Atom::var(GS.eta_left));
        let (right, _) =
            super::CrossSectionGraph::single_th_prefactor_helper_atoms(2, 1, true, true);
        let right_local = right
            .local
            .replace(GS.eta)
            .with(super::Atom::var(GS.eta_right));
        let right_integrated = right
            .integrated
            .replace(GS.eta)
            .with(super::Atom::var(GS.eta_right));

        let legacy_product = ((left_local + left_integrated) * (right_local + right_integrated))
            .expand()
            .replace_multiple(super::CrossSectionGraph::fuse_left_right_replacement());
        let (pieces, legacy) =
            super::CrossSectionGraph::iterated_th_prefactor_helper_atoms(3, 2, 1, 1, true);
        let piece_sum = pieces.local_local
            + pieces.local_integrated
            + pieces.integrated_local
            + pieces.integrated_integrated;

        assert!((legacy_product - &legacy).expand().is_zero());
        assert!((legacy - piece_sum).expand().is_zero());
    }

    #[test]
    fn iterated_threshold_f_params_follow_active_single_axis_shape_order() {
        let params = super::params_for_iterated_threshold_ct(1, 3);
        let f = symbol!("f");
        let base = function!(f, GS.radius_star_left, GS.radius_star_right);
        let expected = [
            base.clone(),
            base.clone().derivative(GS.radius_star_right),
            base.derivative(GS.radius_star_right)
                .derivative(GS.radius_star_right),
        ]
        .into_iter()
        .map(|atom| atom.to_string())
        .collect::<Vec<_>>();
        let actual = params[..expected.len()]
            .iter()
            .map(|atom| atom.to_string())
            .collect::<Vec<_>>();

        assert_eq!(actual, expected);
    }
}