gammalooprs/subtraction/

lu_counterterm.rs

use core::f64;
use std::{collections::BTreeMap, path::Path};

use bincode_trait_derive::{Decode, Encode};
use color_eyre::Result;
use eyre::eyre;
use itertools::Itertools;
use linnet::half_edge::involution::{EdgeIndex, EdgeVec, Orientation};
use spenso::algebra::complex::Complex;
use symbolica::domains::{
    dual::{DualNumberStructure, HyperDual},
    float::{Real, RealLike},
};
use tracing::{debug, warn};
use typed_index_collections::TiVec;

use crate::{
    GammaLoopContext,
    cff::{
        CutCFFIndex,
        esurface::{
            Esurface, EsurfaceCollection, EsurfaceID, ExistingEsurfaceId,
            esurface_value_is_strictly_inside,
        },
        expression::OrientationID,
    },
    graph::{Graph, LmbIndex, LoopMomentumBasis},
    integrands::{
        evaluation::EvaluationMetaData,
        process::{
            GenericEvaluator, ParamBuilder, ThresholdParams,
            evaluators::{
                EvaluatorStack, SingleOrAllOrientations, evaluate_evaluator,
                evaluate_evaluator_single,
            },
            param_builder::LUParams,
        },
    },
    momentum::{
        Rotation, SignOrZero, ThreeMomentum,
        sample::{
            ExternalFourMomenta, ExternalIndex, LoopIndex, LoopMomenta, MomentumSample,
            SubspaceData,
        },
        signature::LoopSignature,
    },
    processes::{
        CutGroupId, EvaluatorBuildTimings, IteratedCtCollection, LUCounterTermData,
        LeftThresholdId, RightThresholdId, build_derivative_structure,
    },
    settings::{GlobalSettings, RuntimeSettings, global::FrozenCompilationMode},
    subtraction::{
        RstarTDependenceEvaluator, RstarTDependenceInput, evaluate_integrated_ct_normalisation,
        evaluate_integrated_ct_normalisation_dual, evaluate_uv_damper, evaluate_uv_damper_dual,
        overlap_subspace::{self, OverlapGroup, OverlapInput, OverlapStructure},
    },
    utils::{
        F, FloatLike,
        hyperdual_utils::{
            DualOrNot, dualize_dual_t_to_dual_r_t, extract_coefficient_t_duals,
            extract_t_derivatives, extract_t_derivatives_complex, new_constant,
            shape_from_cut_cff_index, simple_n_deriv_shape, variable_indices_from_cut_cff_index,
        },
        newton_solver::{NewtonIterationResult, RadialRootDiagnostics, RadialRootIdentity},
    },
};

fn zero_dual_or_not_complex<T: FloatLike>(order: usize, zero: &F<T>) -> DualOrNot<Complex<F<T>>> {
    if order == 0 {
        DualOrNot::NonDual(Complex::new_re(zero.clone()))
    } else {
        DualOrNot::Dual(HyperDual::new(simple_n_deriv_shape(order)))
    }
}

fn negate_dual_or_not_complex<T: FloatLike>(
    value: DualOrNot<Complex<F<T>>>,
) -> DualOrNot<Complex<F<T>>> {
    match value {
        DualOrNot::Dual(dual) => DualOrNot::Dual(-dual),
        DualOrNot::NonDual(non_dual) => DualOrNot::NonDual(-non_dual),
    }
}

fn multiply_dual_or_not_complex<T: FloatLike>(
    lhs: DualOrNot<Complex<F<T>>>,
    rhs: &DualOrNot<Complex<F<T>>>,
) -> DualOrNot<Complex<F<T>>> {
    match (lhs, rhs) {
        (DualOrNot::NonDual(lhs), DualOrNot::NonDual(rhs)) => DualOrNot::NonDual(lhs * rhs),
        (DualOrNot::Dual(lhs), DualOrNot::Dual(rhs)) => DualOrNot::Dual(lhs * rhs.clone()),
        (DualOrNot::Dual(lhs), DualOrNot::NonDual(rhs)) => {
            let rhs_dual = new_constant(&lhs, rhs);
            DualOrNot::Dual(lhs * rhs_dual)
        }
        (DualOrNot::NonDual(lhs), DualOrNot::Dual(rhs)) => {
            let lhs_dual = new_constant(rhs, &lhs);
            DualOrNot::Dual(lhs_dual * rhs.clone())
        }
    }
}

fn plain_t_dual_or_scalar_complex<T: FloatLike>(
    value: &DualOrNot<Complex<F<T>>>,
    t_variable: Option<usize>,
) -> DualOrNot<Complex<F<T>>> {
    match value {
        DualOrNot::Dual(dual) => {
            if let Some(t_variable) = t_variable {
                DualOrNot::Dual(extract_zero_threshold_coefficient_t_dual(dual, t_variable))
            } else {
                DualOrNot::NonDual(dual.values[0].clone())
            }
        }
        DualOrNot::NonDual(value) => DualOrNot::NonDual(value.clone()),
    }
}

fn extract_zero_threshold_coefficient_t_dual<
    T: Clone + spenso::algebra::algebraic_traits::RefZero + Default,
>(
    dual: &HyperDual<T>,
    t_variable: usize,
) -> HyperDual<T> {
    let (coefficient_orders, coefficient_duals) = extract_coefficient_t_duals(dual, t_variable);
    let zero_index = coefficient_orders
        .iter()
        .position(|orders| orders.iter().all(|order| *order == 0))
        .expect("Could not find zero-threshold coefficient in mixed dual");

    coefficient_duals[zero_index].clone()
}

fn embed_t_dual_in_target_shape<T: FloatLike>(
    t_dual: &HyperDual<F<T>>,
    target_shape: &HyperDual<F<T>>,
    t_variable: Option<usize>,
) -> HyperDual<F<T>> {
    match t_variable {
        Some(t_variable) => {
            dualize_dual_t_to_dual_r_t(t_dual.clone(), target_shape.clone(), t_variable)
        }
        None => new_constant(target_shape, &t_dual.values[0]),
    }
}

fn activate_threshold_variable_in_target_shape<T: FloatLike>(
    dual: &mut HyperDual<F<T>>,
    variable: Option<usize>,
) {
    let Some(variable) = variable else {
        return;
    };

    let n_variables = dual.get_shape()[0].len();
    let mut derivative_shape = vec![0; n_variables];
    derivative_shape[variable] = 1;

    let derivative_index = dual
        .get_shape()
        .iter()
        .position(|shape| shape == &derivative_shape)
        .expect("Could not find threshold-variable derivative slot in target shape");

    dual.values[derivative_index] = dual.values[0].one();
}

fn embed_real_dual_or_not_in_target_shape<T: FloatLike>(
    value: &DualOrNot<F<T>>,
    target_shape: &Option<HyperDual<F<T>>>,
    t_variable: Option<usize>,
) -> DualOrNot<F<T>> {
    match (value, target_shape) {
        (DualOrNot::Dual(dual), Some(target_shape)) => {
            DualOrNot::Dual(embed_t_dual_in_target_shape(dual, target_shape, t_variable))
        }
        (DualOrNot::NonDual(value), Some(target_shape)) => {
            DualOrNot::Dual(new_constant(target_shape, value))
        }
        _ => value.clone(),
    }
}

fn embedded_dual_loop_momenta_for_cut_cff_index<T: FloatLike>(
    sample: &MomentumSample<T>,
    target_shape: &Option<HyperDual<F<T>>>,
    t_variable: Option<usize>,
) -> Option<LoopMomenta<HyperDual<F<T>>>> {
    match target_shape {
        Some(target_shape) => Some(match sample.sample.dual_loop_moms.as_ref() {
            Some(dual_loop_moms) => dual_loop_moms
                .iter()
                .map(|momentum| {
                    momentum.map_ref(&|component| {
                        embed_t_dual_in_target_shape(component, target_shape, t_variable)
                    })
                })
                .collect(),
            None => sample
                .loop_moms()
                .iter()
                .map(|momentum| {
                    momentum.map_ref(&|component| new_constant(target_shape, component))
                })
                .collect(),
        }),
        None => sample.sample.dual_loop_moms.clone(),
    }
}

fn extend_plain_helper_real_params<T: FloatLike>(
    params: &mut Vec<Complex<F<T>>>,
    value: &DualOrNot<F<T>>,
    t_variable: Option<usize>,
) {
    match value {
        DualOrNot::Dual(dual) => {
            if let Some(t_variable) = t_variable {
                params.extend(
                    extract_t_derivatives(extract_zero_threshold_coefficient_t_dual(
                        dual, t_variable,
                    ))
                    .into_iter()
                    .map(Complex::new_re),
                );
            } else {
                params.push(Complex::new_re(dual.values[0].clone()));
            }
        }
        DualOrNot::NonDual(value) => params.push(Complex::new_re(value.clone())),
    }
}

fn extend_extracted_f_params<T: FloatLike>(
    params: &mut Vec<Complex<F<T>>>,
    value: &DualOrNot<Complex<F<T>>>,
    t_variable: Option<usize>,
) {
    match value {
        DualOrNot::Dual(dual) => {
            if let Some(t_variable) = t_variable {
                let (_, coefficient_duals) = extract_coefficient_t_duals(dual, t_variable);
                for coefficient_dual in coefficient_duals {
                    params.extend(extract_t_derivatives_complex(coefficient_dual));
                }
            } else {
                params.extend(dual.values.iter().cloned());
            }
        }
        DualOrNot::NonDual(value) => params.push(value.clone()),
    }
}

fn extend_extracted_eta_params<T: FloatLike>(
    params: &mut Vec<Complex<F<T>>>,
    value: &DualOrNot<F<T>>,
    t_variable: Option<usize>,
) {
    match value {
        DualOrNot::Dual(dual) => {
            if let Some(t_variable) = t_variable {
                let (_, coefficient_duals) = extract_coefficient_t_duals(dual, t_variable);
                for coefficient_dual in coefficient_duals {
                    params.extend(
                        extract_t_derivatives(coefficient_dual)
                            .into_iter()
                            .map(Complex::new_re),
                    );
                }
            } else {
                params.extend(dual.values.iter().cloned().map(Complex::new_re));
            }
        }
        DualOrNot::NonDual(value) => params.push(Complex::new_re(value.clone())),
    }
}

fn evaluate_threshold_helper_single<
    T: FloatLike + crate::integrands::process::evaluators::GenericEvaluatorFloat,
>(
    helper: &mut GenericEvaluator,
    cut_cff_index: &CutCFFIndex,
    threshold_result: &DualOrNot<Complex<F<T>>>,
    threshold_params: &ThresholdParams<T>,
    evaluation_meta_data: &mut EvaluationMetaData,
    record_primary_timing: bool,
) -> DualOrNot<Complex<F<T>>> {
    let t_variable = variable_indices_from_cut_cff_index(cut_cff_index).lu_cut;
    let mut helper_params = Vec::new();
    extend_extracted_f_params(&mut helper_params, threshold_result, t_variable);
    extend_extracted_eta_params(
        &mut helper_params,
        &threshold_params.esurface_derivative,
        t_variable,
    );
    extend_plain_helper_real_params(&mut helper_params, &threshold_params.radius, t_variable);
    extend_plain_helper_real_params(
        &mut helper_params,
        &threshold_params.radius_star,
        t_variable,
    );
    extend_plain_helper_real_params(
        &mut helper_params,
        &threshold_params.uv_damp_plus,
        t_variable,
    );
    extend_plain_helper_real_params(
        &mut helper_params,
        &threshold_params.uv_damp_minus,
        t_variable,
    );
    extend_plain_helper_real_params(&mut helper_params, &threshold_params.h_function, t_variable);

    debug!(
        "LU threshold helper input (single): cut_cff_index={:?}, t_variable={:?}, param_count={}",
        cut_cff_index,
        t_variable,
        helper_params.len()
    );
    for (idx, value) in helper_params.iter().enumerate() {
        debug!(
            "LU threshold helper param (single): cut_cff_index={:?}, index={}, value={}",
            cut_cff_index, idx, value
        );
    }

    let mut result = evaluate_evaluator(
        helper,
        &helper_params,
        evaluation_meta_data,
        record_primary_timing,
    );

    debug_assert_eq!(result.len(), 1);
    result.pop().unwrap()
}

fn evaluate_threshold_helper_iterated<
    T: FloatLike + crate::integrands::process::evaluators::GenericEvaluatorFloat,
>(
    helper: &mut GenericEvaluator,
    cut_cff_index: &CutCFFIndex,
    threshold_result: &DualOrNot<Complex<F<T>>>,
    left_threshold_params: &ThresholdParams<T>,
    right_threshold_params: &ThresholdParams<T>,
    evaluation_meta_data: &mut EvaluationMetaData,
    record_primary_timing: bool,
) -> DualOrNot<Complex<F<T>>> {
    let t_variable = variable_indices_from_cut_cff_index(cut_cff_index).lu_cut;
    let mut helper_params = Vec::new();
    extend_extracted_f_params(&mut helper_params, threshold_result, t_variable);
    extend_extracted_eta_params(
        &mut helper_params,
        &left_threshold_params.esurface_derivative,
        t_variable,
    );
    extend_extracted_eta_params(
        &mut helper_params,
        &right_threshold_params.esurface_derivative,
        t_variable,
    );

    extend_plain_helper_real_params(
        &mut helper_params,
        &left_threshold_params.radius,
        t_variable,
    );
    extend_plain_helper_real_params(
        &mut helper_params,
        &left_threshold_params.radius_star,
        t_variable,
    );
    extend_plain_helper_real_params(
        &mut helper_params,
        &left_threshold_params.uv_damp_plus,
        t_variable,
    );
    extend_plain_helper_real_params(
        &mut helper_params,
        &left_threshold_params.uv_damp_minus,
        t_variable,
    );
    extend_plain_helper_real_params(
        &mut helper_params,
        &left_threshold_params.h_function,
        t_variable,
    );

    extend_plain_helper_real_params(
        &mut helper_params,
        &right_threshold_params.radius,
        t_variable,
    );
    extend_plain_helper_real_params(
        &mut helper_params,
        &right_threshold_params.radius_star,
        t_variable,
    );
    extend_plain_helper_real_params(
        &mut helper_params,
        &right_threshold_params.uv_damp_plus,
        t_variable,
    );
    extend_plain_helper_real_params(
        &mut helper_params,
        &right_threshold_params.uv_damp_minus,
        t_variable,
    );
    extend_plain_helper_real_params(
        &mut helper_params,
        &right_threshold_params.h_function,
        t_variable,
    );

    debug!(
        "LU threshold helper input (iterated): cut_cff_index={:?}, t_variable={:?}, param_count={}",
        cut_cff_index,
        t_variable,
        helper_params.len()
    );
    for (idx, value) in helper_params.iter().enumerate() {
        debug!(
            "LU threshold helper param (iterated): cut_cff_index={:?}, index={}, value={}",
            cut_cff_index, idx, value
        );
    }

    let mut result = evaluate_evaluator(
        helper,
        &helper_params,
        evaluation_meta_data,
        record_primary_timing,
    );

    debug_assert_eq!(result.len(), 1);
    result.pop().unwrap()
}

fn real_dual_to_complex<T: FloatLike>(dual: HyperDual<F<T>>) -> HyperDual<Complex<F<T>>> {
    let shape = dual.get_shape().iter().map(|v| v.to_vec()).collect();
    let values = dual.values.into_iter().map(Complex::new_re).collect_vec();
    HyperDual::from_values(shape, values)
}

fn dualize_external_momenta<T: FloatLike>(
    shape: &HyperDual<F<T>>,
    external_moms: &ExternalFourMomenta<F<T>>,
) -> ExternalFourMomenta<HyperDual<F<T>>> {
    external_moms
        .iter()
        .map(|momentum| momentum.map_ref(&|component| new_constant(shape, component)))
        .collect()
}

fn dualize_loop_momenta<T: FloatLike>(
    shape: &HyperDual<F<T>>,
    loop_moms: &LoopMomenta<F<T>>,
) -> LoopMomenta<HyperDual<F<T>>> {
    loop_moms
        .iter()
        .map(|momentum| momentum.map_ref(&|component| new_constant(shape, component)))
        .collect()
}

fn dual_shifted_radius<T: FloatLike>(
    shifted_momenta: &LoopMomenta<HyperDual<F<T>>>,
    subspace: &SubspaceData,
) -> HyperDual<F<T>> {
    let zero = new_constant(
        &shifted_momenta[LoopIndex(0)].px,
        &shifted_momenta[LoopIndex(0)].px.values[0].zero(),
    );

    subspace
        .iter_lmb_indices()
        .fold(zero, |acc, loop_index| {
            acc + shifted_momenta[loop_index].norm_squared()
        })
        .sqrt()
}

fn compute_shift_part_from_dual_momenta_in_subspace<T: FloatLike>(
    esurface: &Esurface,
    loop_moms: &LoopMomenta<HyperDual<F<T>>>,
    external_moms: &ExternalFourMomenta<HyperDual<F<T>>>,
    subspace: &SubspaceData,
    all_lmbs: &TiVec<LmbIndex, LoopMomentumBasis>,
    graph: &Graph,
    masses: &EdgeVec<F<T>>,
) -> HyperDual<F<T>> {
    let lmb = subspace.get_lmb(all_lmbs);
    let zero = new_constant(
        &external_moms[ExternalIndex(0)].temporal.value,
        &external_moms[ExternalIndex(0)].temporal.value.values[0].zero(),
    );

    let full_external_shift = esurface
        .external_shift
        .iter()
        .map(|(index, sign)| {
            let external_signature = &lmb.edge_signatures[*index].external;
            let sign = new_constant(&zero, &F::from_f64(*sign as f64));
            let external_energy = external_signature
                .try_apply(&external_moms.raw)
                .map(|momentum| momentum.temporal.value)
                .unwrap_or_else(|| zero.clone());

            sign * external_energy
        })
        .reduce(|acc, value| acc + value)
        .unwrap_or_else(|| zero.clone());

    let spatial_externals = external_moms
        .iter()
        .map(|momentum| momentum.spatial.clone())
        .collect::<TiVec<ExternalIndex, _>>();

    let remaining_shift = subspace
        .does_not_contain(&esurface.energies, graph)
        .map(|index| {
            let signature = &lmb.edge_signatures[index];
            let momentum = signature
                .try_compute_momentum(&loop_moms.0, &spatial_externals.raw)
                .unwrap_or_else(|| unreachable!());
            let mass = &masses[index];
            let lifted_mass = new_constant(&momentum.px, mass);

            (momentum.norm_squared() + lifted_mass.clone() * lifted_mass).sqrt()
        })
        .reduce(|acc, value| acc + value)
        .unwrap_or_else(|| zero.clone());

    full_external_shift + remaining_shift
}

#[allow(clippy::too_many_arguments)]
fn compute_self_and_r_derivative_subspace_dual<T: FloatLike>(
    esurface: &Esurface,
    radius: &HyperDual<F<T>>,
    shifted_unit_loops_in_subspace: &LoopMomenta<HyperDual<F<T>>>,
    center_in_subspace: &LoopMomenta<HyperDual<F<T>>>,
    external_moms: &ExternalFourMomenta<HyperDual<F<T>>>,
    real_mass_vector: &EdgeVec<F<T>>,
    subspace: &SubspaceData,
    all_lmbs: &TiVec<LmbIndex, LoopMomentumBasis>,
    graph: &Graph,
) -> (HyperDual<F<T>>, HyperDual<F<T>>) {
    let spatial_part_of_externals = external_moms
        .iter()
        .map(|momentum| momentum.spatial.clone())
        .collect::<TiVec<ExternalIndex, _>>();

    let loops: LoopMomenta<HyperDual<F<T>>> = shifted_unit_loops_in_subspace
        .iter_enumerated()
        .map(|(loop_index, shifted_unit_momenta)| {
            if subspace.contains_loop_index(loop_index) {
                shifted_unit_momenta * radius + &center_in_subspace[loop_index]
            } else {
                shifted_unit_momenta.clone()
            }
        })
        .collect();

    let shift = compute_shift_part_from_dual_momenta_in_subspace(
        esurface,
        shifted_unit_loops_in_subspace,
        external_moms,
        subspace,
        all_lmbs,
        graph,
        real_mass_vector,
    );

    let lmb = subspace.get_lmb(all_lmbs);
    let zero = new_constant(radius, &radius.values[0].zero());
    let (derivative, energy_sum) = subspace
        .contains(&esurface.energies, graph)
        .map(|index| {
            let signature = &lmb.edge_signatures[index];
            let momentum = signature
                .try_compute_momentum(&loops.0, &spatial_part_of_externals.raw)
                .unwrap_or_else(|| unreachable!());
            let unit_loop_part = compute_loop_part_subspace_dual(
                &signature.internal,
                shifted_unit_loops_in_subspace,
                subspace,
            );
            let mass = &real_mass_vector[index];
            let lifted_mass = new_constant(&momentum.px, mass);
            let energy = (momentum.norm_squared() + lifted_mass.clone() * lifted_mass).sqrt();
            let numerator = momentum * &unit_loop_part;

            (numerator / &energy, energy)
        })
        .fold((zero.clone(), zero), |(der_sum, en_sum), (der, en)| {
            (der_sum + der, en_sum + en)
        });

    (energy_sum + shift, derivative)
}

fn compute_loop_part_subspace_dual<T: FloatLike>(
    loop_signature: &LoopSignature,
    loop_moms: &LoopMomenta<HyperDual<F<T>>>,
    subspace: &SubspaceData,
) -> ThreeMomentum<HyperDual<F<T>>> {
    let projected: LoopSignature = subspace.project_loop_signature(loop_signature).collect();
    let zero = new_constant(
        &loop_moms[LoopIndex(0)].px,
        &loop_moms[LoopIndex(0)].px.values[0].zero(),
    );
    let mut result = ThreeMomentum::new(zero.clone(), zero.clone(), zero);

    for (loop_index, sign) in projected.iter_enumerated() {
        match sign {
            SignOrZero::Zero => {}
            SignOrZero::Plus => result += &loop_moms[loop_index],
            SignOrZero::Minus => result -= &loop_moms[loop_index],
        }
    }

    result
}

#[derive(Clone, Encode, Decode)]
#[trait_decode(trait = GammaLoopContext)]
pub struct LUThresholdHelperEvaluators {
    pub left_thresholds: TiVec<LeftThresholdId, BTreeMap<CutCFFIndex, GenericEvaluator>>,
    pub right_thresholds: TiVec<RightThresholdId, BTreeMap<CutCFFIndex, GenericEvaluator>>,
    pub iterated: IteratedCtCollection<BTreeMap<CutCFFIndex, GenericEvaluator>>,
}

impl LUThresholdHelperEvaluators {
    fn for_each_generic_evaluator_mut(
        &mut self,
        mut f: impl FnMut(&mut GenericEvaluator) -> color_eyre::Result<()>,
    ) -> color_eyre::Result<()> {
        for evaluators in self.left_thresholds.iter_mut() {
            for evaluator in evaluators.values_mut() {
                f(evaluator)?;
            }
        }
        for evaluators in self.right_thresholds.iter_mut() {
            for evaluator in evaluators.values_mut() {
                f(evaluator)?;
            }
        }
        for evaluators in self.iterated.iter_mut() {
            for evaluator in evaluators.values_mut() {
                f(evaluator)?;
            }
        }

        Ok(())
    }
}

#[derive(Clone, Encode, Decode)]
#[trait_decode(trait = GammaLoopContext)]
pub struct LUCounterTermEvaluators {
    pub left_thresholds_evaluator: TiVec<LeftThresholdId, BTreeMap<CutCFFIndex, EvaluatorStack>>,
    pub right_thresholds_evaluator: TiVec<RightThresholdId, BTreeMap<CutCFFIndex, EvaluatorStack>>,
    pub iterated_evaluator: IteratedCtCollection<BTreeMap<CutCFFIndex, EvaluatorStack>>,
    pub threshold_helpers: LUThresholdHelperEvaluators,
    pub residue_from_e_surface_evaluators: Vec<GenericEvaluator>,
}

impl LUCounterTermEvaluators {
    pub(crate) fn generic_compileable_evaluator_count(&self) -> usize {
        let left = self
            .left_thresholds_evaluator
            .iter()
            .flat_map(|evaluators| evaluators.values())
            .map(EvaluatorStack::generic_evaluator_count)
            .sum::<usize>();
        let right = self
            .right_thresholds_evaluator
            .iter()
            .flat_map(|evaluators| evaluators.values())
            .map(EvaluatorStack::generic_evaluator_count)
            .sum::<usize>();
        let iterated = self
            .iterated_evaluator
            .iter()
            .flat_map(|evaluators| evaluators.values())
            .map(EvaluatorStack::generic_evaluator_count)
            .sum::<usize>();

        // Ignore eager-only helper evaluators and pass-two evaluators.
        left + right + iterated
    }

    pub fn from_atoms(
        counterterm_data: &LUCounterTermData,
        max_cut_order: usize,
        threshold_helpers: LUThresholdHelperEvaluators,
        param_builder: &ParamBuilder,
        settings: &GlobalSettings,
        orientations: &TiVec<OrientationID, EdgeVec<Orientation>>,
    ) -> (Self, EvaluatorBuildTimings) {
        let mut timings = EvaluatorBuildTimings::default();
        let left_thresholds_evaluator = counterterm_data
            .left_atoms
            .iter()
            .map(|parametric_integrands| {
                parametric_integrands
                    .integrands
                    .iter()
                    .map(|(i, atom)| {
                        let dual_shape = shape_from_cut_cff_index(i);

                        let (evaluator, evaluator_timings) = EvaluatorStack::new_with_timings(
                            std::slice::from_ref(atom),
                            param_builder,
                            &orientations.raw,
                            dual_shape,
                            &settings.generation.evaluator,
                        )
                        .unwrap();
                        timings += evaluator_timings;
                        (*i, evaluator)
                    })
                    .collect()
            })
            .collect();

        let right_thresholds_evaluator = counterterm_data
            .right_atoms
            .iter()
            .map(|parametric_integrands| {
                parametric_integrands
                    .integrands
                    .iter()
                    .map(|(i, atom)| {
                        let dual_shape = shape_from_cut_cff_index(i);

                        let (evaluator, evaluator_timings) = EvaluatorStack::new_with_timings(
                            std::slice::from_ref(atom),
                            param_builder,
                            &orientations.raw,
                            dual_shape,
                            &settings.generation.evaluator,
                        )
                        .unwrap();
                        timings += evaluator_timings;
                        (*i, evaluator)
                    })
                    .collect()
            })
            .collect();

        let iterated_timings = std::cell::Cell::new(EvaluatorBuildTimings::default());
        let iterated_evaluator = counterterm_data.iterated.map_ref(|parametric_integrands| {
            parametric_integrands
                .integrands
                .iter()
                .map(|(i, atom)| {
                    let dual_shape = shape_from_cut_cff_index(i);

                    let (evaluator, evaluator_timings) = EvaluatorStack::new_with_timings(
                        std::slice::from_ref(atom),
                        param_builder,
                        &orientations.raw,
                        dual_shape,
                        &settings.generation.evaluator,
                    )
                    .unwrap();
                    let mut timings = iterated_timings.get();
                    timings += evaluator_timings;
                    iterated_timings.set(timings);
                    (*i, evaluator)
                })
                .collect()
        });
        timings += iterated_timings.get();

        let symbolica_started = std::time::Instant::now();
        let pass_two_evaluator = (1..=max_cut_order)
            .map(|order| {
                build_derivative_structure(order as u8, -1, &settings.generation.evaluator)
            })
            .collect();

        timings.symbolica_time += symbolica_started.elapsed();

        (
            LUCounterTermEvaluators {
                left_thresholds_evaluator,
                right_thresholds_evaluator,
                iterated_evaluator,
                threshold_helpers,
                residue_from_e_surface_evaluators: pass_two_evaluator,
            },
            timings,
        )
    }

    pub(crate) fn compile(
        &mut self,
        path: impl AsRef<Path>,
        cut_group_id: CutGroupId,
        frozen_mode: &FrozenCompilationMode,
    ) -> color_eyre::Result<()> {
        for (threshold_id, evaluators) in self.left_thresholds_evaluator.iter_mut_enumerated() {
            for (index, evaluator) in evaluators.iter_mut() {
                let name = format!(
                    "cut_group_{}_left_threshold_{}_index_{}",
                    cut_group_id.0, threshold_id.0, index
                );
                evaluator.compile(&name, path.as_ref(), frozen_mode)?;
            }
        }

        for (threshold_id, evaluators) in self.right_thresholds_evaluator.iter_mut_enumerated() {
            for (index, evaluator) in evaluators.iter_mut() {
                let name = format!(
                    "cut_group_{}_right_threshold_{}_index_{}",
                    cut_group_id.0, threshold_id.0, index
                );
                evaluator.compile(&name, path.as_ref(), frozen_mode)?;
            }
        }

        for (iterated_index, evaluators) in self.iterated_evaluator.iter_mut().enumerate() {
            for (index, evaluator) in evaluators.iter_mut() {
                let name = format!(
                    "cut_group_{}_iterated_{}_index_{}",
                    cut_group_id.0, iterated_index, index
                );
                evaluator.compile(&name, path.as_ref(), frozen_mode)?;
            }
        }

        for (order, pass_to_evaluator) in self
            .residue_from_e_surface_evaluators
            .iter_mut()
            .enumerate()
        {
            let name = format!("cut_group_{}_pass_two_{}", cut_group_id.0, order);
            pass_to_evaluator.compile_external(
                path.as_ref().join(&name).with_extension("cpp"),
                &name,
                path.as_ref().join(&name).with_extension("so"),
                frozen_mode,
            )?;
        }

        Ok(())
    }

    pub(crate) fn for_each_generic_evaluator_mut(
        &mut self,
        mut f: impl FnMut(&mut GenericEvaluator) -> color_eyre::Result<()>,
    ) -> color_eyre::Result<()> {
        for evaluators in self.left_thresholds_evaluator.iter_mut() {
            for evaluator in evaluators.values_mut() {
                evaluator.for_each_generic_evaluator_mut(&mut f)?;
            }
        }
        for evaluators in self.right_thresholds_evaluator.iter_mut() {
            for evaluator in evaluators.values_mut() {
                evaluator.for_each_generic_evaluator_mut(&mut f)?;
            }
        }
        for evaluators in self.iterated_evaluator.iter_mut() {
            for evaluator in evaluators.values_mut() {
                evaluator.for_each_generic_evaluator_mut(&mut f)?;
            }
        }
        self.threshold_helpers
            .for_each_generic_evaluator_mut(&mut f)?;
        for pass_to_evaluator in self.residue_from_e_surface_evaluators.iter_mut() {
            f(pass_to_evaluator)?;
        }

        Ok(())
    }
}

type CutThresholds = (
    TiVec<LeftThresholdId, Esurface>,
    TiVec<RightThresholdId, Esurface>,
);

#[derive(Clone, Encode, Decode)]
#[trait_decode(trait = GammaLoopContext)]
pub(crate) struct LUCounterTerm {
    pub evaluators: TiVec<CutGroupId, LUCounterTermEvaluators>,
    pub thresholds: TiVec<CutGroupId, CutThresholds>,
    pub subspaces: TiVec<CutGroupId, (SubspaceData, SubspaceData)>,
    pub active_cut_groups: TiVec<CutGroupId, bool>,
    pub active_left_thresholds: TiVec<CutGroupId, TiVec<LeftThresholdId, bool>>,
    pub active_right_thresholds: TiVec<CutGroupId, TiVec<RightThresholdId, bool>>,
    pub active_iterated_thresholds: TiVec<CutGroupId, IteratedCtCollection<bool>>,
    pub rstar_dependence_calculator: TiVec<CutGroupId, RstarTDependenceEvaluator>,
}

pub struct LUCTKinematicPoint<T: FloatLike> {
    pub unrescaled_sample: MomentumSample<T>,
    pub dualized_momentum_sample_cache: Vec<MomentumSample<T>>,
    pub lu_cut_parameter_cache: Vec<LUParams<T>>,
    pub lu_cut_esurface_values: Vec<DualOrNot<F<T>>>,
}

impl<T: FloatLike> LUCTKinematicPoint<T> {
    pub fn unrescaled_sample(&self) -> &MomentumSample<T> {
        &self.unrescaled_sample
    }

    pub fn representative_sample(&self) -> &MomentumSample<T> {
        &self.dualized_momentum_sample_cache[0]
    }

    pub fn sample_for_order(&self, order: usize) -> &MomentumSample<T> {
        &self.dualized_momentum_sample_cache[order]
    }

    pub fn lmb_transform(&self, from: &LoopMomentumBasis, to: &LoopMomentumBasis) -> Self {
        let transformed_samples = self
            .dualized_momentum_sample_cache
            .iter()
            .map(|sample| sample.lmb_transform(from, to))
            .collect();

        Self {
            unrescaled_sample: self.unrescaled_sample.clone(),
            dualized_momentum_sample_cache: transformed_samples,
            lu_cut_parameter_cache: self.lu_cut_parameter_cache.clone(),
            lu_cut_esurface_values: self.lu_cut_esurface_values.clone(),
        }
    }

    pub fn non_dual_cut_params(&self) -> LUParams<T> {
        self.lu_cut_parameter_cache[0].clone()
    }

    pub fn cut_params_for_order(&self, order: usize) -> &LUParams<T> {
        &self.lu_cut_parameter_cache[order]
    }

    pub fn cut_params_for_cut_cff_index(&self, cut_cff_index: &CutCFFIndex) -> LUParams<T> {
        let lu_order = cut_cff_index
            .lu_cut_order
            .expect("LU cut parameter lookup requires lu_cut_order")
            - 1;
        let base_params = self.cut_params_for_order(lu_order).clone();
        let target_shape = shape_from_cut_cff_index(cut_cff_index).map(HyperDual::new);
        let t_variable = variable_indices_from_cut_cff_index(cut_cff_index).lu_cut;

        LUParams {
            tstar: embed_real_dual_or_not_in_target_shape(
                &base_params.tstar,
                &target_shape,
                t_variable,
            ),
            h_function: embed_real_dual_or_not_in_target_shape(
                &base_params.h_function,
                &target_shape,
                t_variable,
            ),
        }
    }

    pub fn cut_esurface_for_order(&self, order: usize) -> &DualOrNot<F<T>> {
        &self.lu_cut_esurface_values[order]
    }

    pub fn non_dual_cut_esurface_value(&self) -> F<T> {
        match self.cut_esurface_for_order(0) {
            DualOrNot::NonDual(value) => value.clone(),
            DualOrNot::Dual(_) => {
                unreachable!("representative LU cut e-surface value must be non-dual")
            }
        }
    }

    pub fn new(unrescaled_sample: MomentumSample<T>) -> Self {
        Self {
            unrescaled_sample,
            dualized_momentum_sample_cache: Vec::new(),
            lu_cut_parameter_cache: Vec::new(),
            lu_cut_esurface_values: Vec::new(),
        }
    }
}

impl LUCounterTerm {
    fn radial_root_identity(
        graph_name: &str,
        cut_group_id: CutGroupId,
        side: &str,
        overlap_group: usize,
        esurface_id: EsurfaceID,
        probe_rotation: &Rotation,
    ) -> RadialRootIdentity {
        RadialRootIdentity::new(format!(
            "LU graph '{graph_name}' cut group {} {side} overlap group {overlap_group} E-surface {} probe rotation {}",
            cut_group_id.0, esurface_id.0, probe_rotation.method,
        ))
    }

    pub(crate) fn cut_group_is_active(&self, cut_group_id: CutGroupId) -> bool {
        self.active_cut_groups[cut_group_id]
    }

    pub(crate) fn compile(
        &mut self,
        path: impl AsRef<Path>,
        frozen_mode: &FrozenCompilationMode,
    ) -> color_eyre::Result<()> {
        for (cut_group_id, evaluators) in self.evaluators.iter_mut_enumerated() {
            evaluators.compile(path.as_ref(), cut_group_id, frozen_mode)?;
        }
        Ok(())
    }

    pub(crate) fn for_each_generic_evaluator_mut(
        &mut self,
        mut f: impl FnMut(&mut GenericEvaluator) -> color_eyre::Result<()>,
    ) -> color_eyre::Result<()> {
        for evaluators in self.evaluators.iter_mut() {
            evaluators.for_each_generic_evaluator_mut(&mut f)?;
        }
        Ok(())
    }

    fn ensure_active_cut_group(&self, cut_group_id: CutGroupId) -> Result<()> {
        if self.cut_group_is_active(cut_group_id) {
            return Ok(());
        }

        Err(eyre!(
            "Cut group {} was reached at runtime even though generation marked it inactive for the selected orientation subset",
            cut_group_id.0
        ))
    }

    fn ensure_active_left_threshold(
        &self,
        cut_group_id: CutGroupId,
        threshold_id: LeftThresholdId,
    ) -> Result<()> {
        if self.active_left_thresholds[cut_group_id][threshold_id] {
            return Ok(());
        }

        Err(eyre!(
            "Left threshold evaluator {} for cut group {} was reached at runtime even though generation marked it inactive for the selected orientation subset",
            threshold_id.0,
            cut_group_id.0
        ))
    }

    fn ensure_active_right_threshold(
        &self,
        cut_group_id: CutGroupId,
        threshold_id: RightThresholdId,
    ) -> Result<()> {
        if self.active_right_thresholds[cut_group_id][threshold_id] {
            return Ok(());
        }

        Err(eyre!(
            "Right threshold evaluator {} for cut group {} was reached at runtime even though generation marked it inactive for the selected orientation subset",
            threshold_id.0,
            cut_group_id.0
        ))
    }

    fn ensure_active_iterated_threshold(
        &self,
        cut_group_id: CutGroupId,
        iterated_index: (LeftThresholdId, RightThresholdId),
    ) -> Result<()> {
        if self.active_iterated_thresholds[cut_group_id][iterated_index] {
            return Ok(());
        }

        Err(eyre!(
            "Iterated threshold evaluator ({}, {}) for cut group {} was reached at runtime even though generation marked it inactive for the selected orientation subset",
            iterated_index.0.0,
            iterated_index.1.0,
            cut_group_id.0
        ))
    }

    #[allow(clippy::too_many_arguments)]
    pub(crate) fn evaluate<T: FloatLike>(
        &mut self,
        kinematic_point: &LUCTKinematicPoint<T>,
        cut_group_id: CutGroupId,
        reversed_edges: &[EdgeIndex],
        all_lmbs: &TiVec<LmbIndex, LoopMomentumBasis>,
        graph: &Graph,
        masses: &EdgeVec<F<T>>,
        probe_rotation: &Rotation,
        settings: &RuntimeSettings,
        param_builder: &mut ParamBuilder<f64>,
        orientations: SingleOrAllOrientations<'_, OrientationID>,
        evaluation_meta_data: &mut EvaluationMetaData,
        record_primary_timing: bool,
    ) -> Result<Complex<F<T>>> {
        self.ensure_active_cut_group(cut_group_id)?;

        let (left_thresholds_typed, right_thresholds_typed) = &self.thresholds[cut_group_id];

        let left_thresholds = TiVec::from_ref(&left_thresholds_typed.raw);
        let right_thresholds = TiVec::from_ref(&right_thresholds_typed.raw);

        if left_thresholds.is_empty() && right_thresholds.is_empty() {
            return Ok(Complex::new_re(
                kinematic_point.representative_sample().zero(),
            ));
        }

        let (left_subspace, right_subspace) = &self.subspaces[cut_group_id];
        let (sample_left_transformed, sample_right_transformed) = (
            kinematic_point
                .lmb_transform(&graph.loop_momentum_basis, left_subspace.get_lmb(all_lmbs)),
            kinematic_point
                .lmb_transform(&graph.loop_momentum_basis, right_subspace.get_lmb(all_lmbs)),
        );

        debug!("possible left thresholds: {}", left_thresholds.len());
        debug!("possible right thresholds: {}", right_thresholds.len());

        let masses_f64: EdgeVec<F<f64>> = masses.iter().map(|(_, m)| F(m.to_f64())).collect();
        let sample_left_transformed_f64 = sample_left_transformed
            .representative_sample()
            .loop_moms()
            .iter()
            .map(|lm| lm.to_f64())
            .collect();
        let sample_right_transformed_f64 = sample_right_transformed
            .representative_sample()
            .loop_moms()
            .iter()
            .map(|lm| lm.to_f64())
            .collect();
        let external_moms_f64 = kinematic_point
            .representative_sample()
            .external_moms()
            .iter()
            .map(|em| em.to_f64())
            .collect();

        let e_cm = F::from_f64(settings.kinematics.e_cm);
        let esurface_existence_threshold =
            F::from_f64(settings.subtraction.esurface_existence_threshold);

        let left_existing_esurfaces = self.thresholds[cut_group_id]
            .0
            .iter_enumerated()
            .filter_map(|(left_id, esurface)| {
                let classification = esurface.classify_existence_subspace(
                    sample_left_transformed.representative_sample().loop_moms(),
                    kinematic_point.representative_sample().external_moms(),
                    left_subspace,
                    all_lmbs,
                    graph,
                    masses,
                    reversed_edges,
                    &e_cm,
                    &esurface_existence_threshold,
                );
                debug!(
                    graph = %graph.name,
                    side = "left",
                    esurface_id = left_id.0,
                    status = classification.label(),
                    normalized_margin = ?classification.normalized_margin(),
                    non_existing_reason = ?classification.non_existing_reason(),
                    "classified LU threshold surface"
                );
                if classification.is_existing() {
                    Some(EsurfaceID::from(left_id.0))
                } else {
                    None
                }
            })
            .collect::<TiVec<ExistingEsurfaceId, _>>();

        let right_existing_esurfaces = self.thresholds[cut_group_id]
            .1
            .iter_enumerated()
            .filter_map(|(right_id, esurface)| {
                let classification = esurface.classify_existence_subspace(
                    sample_right_transformed.representative_sample().loop_moms(),
                    kinematic_point.representative_sample().external_moms(),
                    right_subspace,
                    all_lmbs,
                    graph,
                    masses,
                    reversed_edges,
                    &e_cm,
                    &esurface_existence_threshold,
                );
                debug!(
                    graph = %graph.name,
                    side = "right",
                    esurface_id = right_id.0,
                    status = classification.label(),
                    normalized_margin = ?classification.normalized_margin(),
                    non_existing_reason = ?classification.non_existing_reason(),
                    "classified LU threshold surface"
                );
                if classification.is_existing() {
                    Some(EsurfaceID::from(right_id.0))
                } else {
                    None
                }
            })
            .collect::<TiVec<ExistingEsurfaceId, _>>();

        debug!(
            "number of thresholds on the left: {}",
            left_existing_esurfaces.len()
        );
        debug!(
            "number of thresholds on the right: {}",
            right_existing_esurfaces.len()
        );

        if left_existing_esurfaces.is_empty() && right_existing_esurfaces.is_empty() {
            return Ok(Complex::new_re(
                kinematic_point.representative_sample().zero(),
            ));
        }

        let left_overlap_input = OverlapInput {
            graph,
            subspace: left_subspace,
            settings,
            edge_masses: masses_f64.clone(),
            lmbs: all_lmbs,
            thresholds: left_thresholds,
        };

        let right_overlap_input = OverlapInput {
            graph,
            subspace: right_subspace,
            settings,
            edge_masses: masses_f64,
            lmbs: all_lmbs,
            thresholds: right_thresholds,
        };

        let left_overlap = match overlap_subspace::find_maximal_overlap(
            &left_overlap_input,
            &left_existing_esurfaces,
            &sample_left_transformed_f64,
            &external_moms_f64,
            probe_rotation,
        ) {
            Ok(left_overlap) => left_overlap,
            Err(error) => {
                evaluation_meta_data.record_threshold_counterterm_error(format!(
                    "LU graph '{}' cut group {} left subspace failed overlap-center construction in probe rotation {}: {}",
                    graph.name,
                    cut_group_id.0,
                    probe_rotation.method,
                    error,
                ));
                return Ok(Complex::new_re(F::from_f64(f64::NAN)));
            }
        };

        let right_overlap = match overlap_subspace::find_maximal_overlap(
            &right_overlap_input,
            &right_existing_esurfaces,
            &sample_right_transformed_f64,
            &external_moms_f64,
            probe_rotation,
        ) {
            Ok(right_overlap) => right_overlap,
            Err(error) => {
                evaluation_meta_data.record_threshold_counterterm_error(format!(
                    "LU graph '{}' cut group {} right subspace failed overlap-center construction in probe rotation {}: {}",
                    graph.name,
                    cut_group_id.0,
                    probe_rotation.method,
                    error,
                ));
                return Ok(Complex::new_re(F::from_f64(f64::NAN)));
            }
        };

        debug!("left overlap structure: {}", left_overlap);

        debug!("right overlap structure: {}", right_overlap);

        let left_counterterm_builder = CounterTermBuilder::new(
            graph,
            settings,
            left_overlap_input.thresholds,
            sample_left_transformed,
            &left_overlap,
            masses,
            all_lmbs,
            left_subspace,
            probe_rotation,
        );

        let left_overlap_builders = left_overlap
            .overlap_groups
            .iter()
            .map(|overlap_group| left_counterterm_builder.new_overlap_builder(overlap_group))
            .collect_vec();

        let mut radial_root_failed = false;
        let mut left_overlap_solutions = Vec::with_capacity(left_overlap_builders.len());
        for (overlap_group_index, overlap_builder) in left_overlap_builders.iter().enumerate() {
            let mut group_solutions =
                Vec::with_capacity(overlap_builder.overlap_group.existing_esurfaces.len());
            for &existing_esurface_id in &overlap_builder.overlap_group.existing_esurfaces {
                let esurface_id = left_overlap.existing_esurfaces[existing_esurface_id];
                let radial_root_identity = Self::radial_root_identity(
                    &graph.name,
                    cut_group_id,
                    "left",
                    overlap_group_index,
                    esurface_id,
                    probe_rotation,
                );
                let Some(solution) = overlap_builder
                    .new_esurface_builder(existing_esurface_id)
                    .solve_rstar(
                        &mut self.rstar_dependence_calculator[cut_group_id],
                        &radial_root_identity,
                        &mut evaluation_meta_data.radial_root_diagnostics,
                    )
                else {
                    evaluation_meta_data.record_threshold_counterterm_error(format!(
                        "LU graph '{}' cut group {} left overlap group {} E-surface {} failed center or radial-root validation in probe rotation {}",
                        graph.name,
                        cut_group_id.0,
                        overlap_group_index,
                        esurface_id.0,
                        probe_rotation.method,
                    ));
                    radial_root_failed = true;
                    continue;
                };
                group_solutions.push(solution);
            }
            left_overlap_solutions.push(group_solutions);
        }

        let right_counterterm_builder = CounterTermBuilder::new(
            graph,
            settings,
            right_overlap_input.thresholds,
            sample_right_transformed,
            &right_overlap,
            masses,
            all_lmbs,
            right_subspace,
            probe_rotation,
        );

        let right_overlap_builders = right_overlap
            .overlap_groups
            .iter()
            .map(|overlap_group| right_counterterm_builder.new_overlap_builder(overlap_group))
            .collect_vec();

        let mut right_overlap_solutions = Vec::with_capacity(right_overlap_builders.len());
        for (overlap_group_index, overlap_builder) in right_overlap_builders.iter().enumerate() {
            let mut group_solutions =
                Vec::with_capacity(overlap_builder.overlap_group.existing_esurfaces.len());
            for &existing_esurface_id in &overlap_builder.overlap_group.existing_esurfaces {
                let esurface_id = right_overlap.existing_esurfaces[existing_esurface_id];
                let radial_root_identity = Self::radial_root_identity(
                    &graph.name,
                    cut_group_id,
                    "right",
                    overlap_group_index,
                    esurface_id,
                    probe_rotation,
                );
                let Some(solution) = overlap_builder
                    .new_esurface_builder(existing_esurface_id)
                    .solve_rstar(
                        &mut self.rstar_dependence_calculator[cut_group_id],
                        &radial_root_identity,
                        &mut evaluation_meta_data.radial_root_diagnostics,
                    )
                else {
                    evaluation_meta_data.record_threshold_counterterm_error(format!(
                        "LU graph '{}' cut group {} right overlap group {} E-surface {} failed center or radial-root validation in probe rotation {}",
                        graph.name,
                        cut_group_id.0,
                        overlap_group_index,
                        esurface_id.0,
                        probe_rotation.method,
                    ));
                    radial_root_failed = true;
                    continue;
                };
                group_solutions.push(solution);
            }
            right_overlap_solutions.push(group_solutions);
        }

        if radial_root_failed {
            return Ok(Complex::new_re(F::from_f64(f64::NAN)));
        }

        let zero = kinematic_point.representative_sample().zero();
        let mut total_result = Complex::new_re(zero.clone());

        for order in 0..kinematic_point.lu_cut_parameter_cache.len() {
            let mut left_evaluations = zero_dual_or_not_complex(order, &zero);

            for solutions_group in left_overlap_solutions.iter() {
                for solution in solutions_group {
                    let representative_sample = solution.rstar_sample_for_order(order);
                    let left_threshold_id =
                        LeftThresholdId::from(representative_sample.get_esurface_id().0);
                    self.ensure_active_left_threshold(cut_group_id, left_threshold_id)?;

                    let matching_cut_indices = self.evaluators[cut_group_id]
                        .left_thresholds_evaluator[left_threshold_id]
                        .keys()
                        .filter(|cut_cff_index| {
                            cut_cff_index.lu_cut_order == Some(order + 1)
                                && cut_cff_index.right_threshold_order.is_none()
                        })
                        .copied()
                        .collect_vec();

                    for cut_cff_index in matching_cut_indices {
                        let variable_indices = variable_indices_from_cut_cff_index(&cut_cff_index);
                        let lu_cut_params =
                            kinematic_point.cut_params_for_cut_cff_index(&cut_cff_index);
                        let sample = solution.rstar_sample_for_cut_cff_index(
                            &cut_cff_index,
                            variable_indices.left_threshold,
                        );
                        debug!("left threshold parameters");

                        let left_threshold_params: ThresholdParams<T> =
                            sample.extract_threshold_parameters(true);
                        debug!(
                            "LU left evaluator input: cut_group_id={}, left_threshold_id={}, cut_cff_index={:?}, lu_order={}, r={}, rstar={}",
                            cut_group_id.0,
                            left_threshold_id.0,
                            cut_cff_index,
                            order + 1,
                            left_threshold_params.radius,
                            left_threshold_params.radius_star,
                        );
                        let inverse_transformed_sample = sample.get_inverse_transformed_sample();

                        let params = T::get_parameters(
                            param_builder,
                            (false, false),
                            graph,
                            &inverse_transformed_sample,
                            settings.kinematics.externals.get_helicities(),
                            &settings.additional_params(),
                            Some(&left_threshold_params),
                            None,
                            Some(&lu_cut_params),
                        );

                        let result_of_this_ct = self.evaluators[cut_group_id]
                            .left_thresholds_evaluator[left_threshold_id]
                            .get_mut(&cut_cff_index)
                            .unwrap()
                            .evaluate(
                                params,
                                orientations,
                                settings,
                                evaluation_meta_data,
                                record_primary_timing,
                            )
                            .unwrap()
                            .pop()
                            .unwrap();

                        debug!(
                            "result of left threshold evaluator {:?}: {}",
                            cut_cff_index, result_of_this_ct
                        );

                        let helper_completed_result = evaluate_threshold_helper_single(
                            self.evaluators[cut_group_id]
                                .threshold_helpers
                                .left_thresholds[left_threshold_id]
                                .get_mut(&cut_cff_index)
                                .unwrap(),
                            &cut_cff_index,
                            &result_of_this_ct,
                            &left_threshold_params,
                            evaluation_meta_data,
                            record_primary_timing,
                        );

                        let multi_channeling_factor = plain_t_dual_or_scalar_complex(
                            &sample.value_of_multi_channeling_factor,
                            variable_indices.lu_cut,
                        );

                        left_evaluations += multiply_dual_or_not_complex(
                            helper_completed_result,
                            &multi_channeling_factor,
                        );
                    }
                }
            }

            let mut right_evaluations = zero_dual_or_not_complex(order, &zero);

            for solutions_group in right_overlap_solutions.iter() {
                for solution in solutions_group {
                    let representative_sample = solution.rstar_sample_for_order(order);
                    let right_threshold_id =
                        RightThresholdId::from(representative_sample.get_esurface_id().0);
                    self.ensure_active_right_threshold(cut_group_id, right_threshold_id)?;

                    let matching_cut_indices = self.evaluators[cut_group_id]
                        .right_thresholds_evaluator[right_threshold_id]
                        .keys()
                        .filter(|cut_cff_index| {
                            cut_cff_index.lu_cut_order == Some(order + 1)
                                && cut_cff_index.left_threshold_order.is_none()
                        })
                        .copied()
                        .collect_vec();

                    for cut_cff_index in matching_cut_indices {
                        let variable_indices = variable_indices_from_cut_cff_index(&cut_cff_index);
                        let lu_cut_params =
                            kinematic_point.cut_params_for_cut_cff_index(&cut_cff_index);
                        let sample = solution.rstar_sample_for_cut_cff_index(
                            &cut_cff_index,
                            variable_indices.right_threshold,
                        );
                        debug!("right threshold parameters");
                        let right_threshold_params: ThresholdParams<T> =
                            sample.extract_threshold_parameters(true);
                        debug!(
                            "LU right evaluator input: cut_group_id={}, right_threshold_id={}, cut_cff_index={:?}, lu_order={}, r={}, rstar={}",
                            cut_group_id.0,
                            right_threshold_id.0,
                            cut_cff_index,
                            order + 1,
                            right_threshold_params.radius,
                            right_threshold_params.radius_star,
                        );
                        let inverse_transformed_sample = sample.get_inverse_transformed_sample();

                        let params = T::get_parameters(
                            param_builder,
                            (false, false),
                            graph,
                            &inverse_transformed_sample,
                            settings.kinematics.externals.get_helicities(),
                            &settings.additional_params(),
                            None,
                            Some(&right_threshold_params),
                            Some(&lu_cut_params),
                        );

                        let result_of_this_ct = self.evaluators[cut_group_id]
                            .right_thresholds_evaluator[right_threshold_id]
                            .get_mut(&cut_cff_index)
                            .unwrap()
                            .evaluate(
                                params,
                                orientations,
                                settings,
                                evaluation_meta_data,
                                record_primary_timing,
                            )
                            .unwrap()
                            .pop()
                            .unwrap();

                        debug!(
                            "result of right threshold evaluator {:?}: {}",
                            cut_cff_index, result_of_this_ct
                        );

                        let helper_completed_result = evaluate_threshold_helper_single(
                            self.evaluators[cut_group_id]
                                .threshold_helpers
                                .right_thresholds[right_threshold_id]
                                .get_mut(&cut_cff_index)
                                .unwrap(),
                            &cut_cff_index,
                            &result_of_this_ct,
                            &right_threshold_params,
                            evaluation_meta_data,
                            record_primary_timing,
                        );

                        let multi_channeling_factor = plain_t_dual_or_scalar_complex(
                            &sample.value_of_multi_channeling_factor,
                            variable_indices.lu_cut,
                        );

                        right_evaluations += multiply_dual_or_not_complex(
                            helper_completed_result,
                            &multi_channeling_factor,
                        );
                    }
                }
            }

            let mut cartesian_product_result = zero_dual_or_not_complex(order, &zero);

            for left_solutions_group in left_overlap_solutions.iter() {
                for left_solution in left_solutions_group {
                    let left_sample_representative = left_solution.rstar_sample_for_order(order);
                    for right_solutions_group in right_overlap_solutions.iter() {
                        for right_solution in right_solutions_group {
                            let right_sample_representative =
                                right_solution.rstar_sample_for_order(order);
                            let iterated_index = (
                                LeftThresholdId::from(
                                    left_sample_representative.get_esurface_id().0,
                                ),
                                RightThresholdId::from(
                                    right_sample_representative.get_esurface_id().0,
                                ),
                            );
                            self.ensure_active_iterated_threshold(cut_group_id, iterated_index)?;

                            let matching_cut_indices = self.evaluators[cut_group_id]
                                .iterated_evaluator[iterated_index]
                                .keys()
                                .filter(|cut_cff_index| {
                                    cut_cff_index.lu_cut_order == Some(order + 1)
                                })
                                .copied()
                                .collect_vec();

                            for cut_cff_index in matching_cut_indices {
                                let variable_indices =
                                    variable_indices_from_cut_cff_index(&cut_cff_index);
                                let lu_cut_params =
                                    kinematic_point.cut_params_for_cut_cff_index(&cut_cff_index);
                                let sample_left = left_solution.rstar_sample_for_cut_cff_index(
                                    &cut_cff_index,
                                    variable_indices.left_threshold,
                                );
                                let sample_right = right_solution.rstar_sample_for_cut_cff_index(
                                    &cut_cff_index,
                                    variable_indices.right_threshold,
                                );

                                let left_threshold_params: ThresholdParams<T> =
                                    sample_left.extract_threshold_parameters(false);
                                let right_threshold_params: ThresholdParams<T> =
                                    sample_right.extract_threshold_parameters(false);
                                debug!(
                                    "LU iterated evaluator input: cut_group_id={}, left_threshold_id={}, right_threshold_id={}, cut_cff_index={:?}, lu_order={}, left_r={}, left_rstar={}, right_r={}, right_rstar={}",
                                    cut_group_id.0,
                                    iterated_index.0.0,
                                    iterated_index.1.0,
                                    cut_cff_index,
                                    order + 1,
                                    left_threshold_params.radius,
                                    left_threshold_params.radius_star,
                                    right_threshold_params.radius,
                                    right_threshold_params.radius_star,
                                );
                                let multi_channeling_factor = plain_t_dual_or_scalar_complex(
                                    &multiply_dual_or_not_complex(
                                        sample_left.value_of_multi_channeling_factor.clone(),
                                        &sample_right.value_of_multi_channeling_factor,
                                    ),
                                    variable_indices.lu_cut,
                                );
                                let inverse_transformed_momentum_sample =
                                    merge_and_inverse_transform(&sample_left, &sample_right);

                                let params = T::get_parameters(
                                    param_builder,
                                    (false, false),
                                    graph,
                                    &inverse_transformed_momentum_sample,
                                    settings.kinematics.externals.get_helicities(),
                                    &settings.additional_params(),
                                    Some(&left_threshold_params),
                                    Some(&right_threshold_params),
                                    Some(&lu_cut_params),
                                );

                                let result_of_this_ct = self.evaluators[cut_group_id]
                                    .iterated_evaluator[iterated_index]
                                    .get_mut(&cut_cff_index)
                                    .unwrap()
                                    .evaluate(
                                        params,
                                        orientations,
                                        settings,
                                        evaluation_meta_data,
                                        record_primary_timing,
                                    )
                                    .unwrap()
                                    .pop()
                                    .unwrap();

                                let helper_completed_result = evaluate_threshold_helper_iterated(
                                    self.evaluators[cut_group_id].threshold_helpers.iterated
                                        [iterated_index]
                                        .get_mut(&cut_cff_index)
                                        .unwrap(),
                                    &cut_cff_index,
                                    &result_of_this_ct,
                                    &left_threshold_params,
                                    &right_threshold_params,
                                    evaluation_meta_data,
                                    record_primary_timing,
                                );

                                cartesian_product_result += multiply_dual_or_not_complex(
                                    helper_completed_result,
                                    &multi_channeling_factor,
                                );
                            }
                        }
                    }
                }
            }

            let mut pass_one_result = left_evaluations.clone();
            pass_one_result += right_evaluations;
            pass_one_result += negate_dual_or_not_complex(cartesian_product_result);
            let pass_one_result = negate_dual_or_not_complex(pass_one_result);

            let mut params_for_pass_two = vec![];
            match pass_one_result {
                DualOrNot::Dual(dual_result) => {
                    params_for_pass_two
                        .extend_from_slice(&extract_t_derivatives_complex(dual_result));
                }
                DualOrNot::NonDual(non_dual_result) => {
                    params_for_pass_two.push(non_dual_result);
                }
            }

            match kinematic_point.cut_esurface_for_order(order).clone() {
                DualOrNot::Dual(dual_e_surface) => {
                    extract_t_derivatives(dual_e_surface)[1..]
                        .iter()
                        .for_each(|value| params_for_pass_two.push(Complex::new_re(value.clone())));
                }
                DualOrNot::NonDual(non_dual_e_surface) => {
                    debug!("non dual esurface: {}", non_dual_e_surface);
                    params_for_pass_two.push(Complex::new_re(non_dual_e_surface));
                }
            }

            debug!(
                "LU pass-two evaluator input: cut_group_id={}, lu_order={}, param_count={}",
                cut_group_id.0,
                order + 1,
                params_for_pass_two.len()
            );
            for (idx, value) in params_for_pass_two.iter().enumerate() {
                debug!(
                    "LU pass-two evaluator param: cut_group_id={}, lu_order={}, index={}, value={}",
                    cut_group_id.0,
                    order + 1,
                    idx,
                    value
                );
            }

            let pass_two_result = evaluate_evaluator_single(
                &mut self.evaluators[cut_group_id].residue_from_e_surface_evaluators[order],
                &params_for_pass_two,
                evaluation_meta_data,
                record_primary_timing,
            );

            total_result += pass_two_result;
        }

        Ok(total_result)
    }
}

struct CounterTermBuilder<'a, T: FloatLike> {
    overlap_structure: &'a OverlapStructure,
    real_mass_vector: &'a EdgeVec<F<T>>,
    e_cm: F<T>,
    graph: &'a Graph,
    subspace: &'a SubspaceData,
    all_lmbs: &'a TiVec<LmbIndex, LoopMomentumBasis>,
    settings: &'a RuntimeSettings,
    esurface_collection: &'a EsurfaceCollection,
    transformed_kinematic_point: LUCTKinematicPoint<T>,
    probe_rotation: &'a Rotation,
}

impl<'a, T: FloatLike> CounterTermBuilder<'a, T> {
    #[allow(clippy::too_many_arguments)]
    fn new(
        graph: &'a Graph,
        settings: &'a RuntimeSettings,
        esurface_collection: &'a EsurfaceCollection,
        transformed_kinematic_point: LUCTKinematicPoint<T>,
        overlap_structure: &'a OverlapStructure,
        masses: &'a EdgeVec<F<T>>,
        all_lmbs: &'a TiVec<LmbIndex, LoopMomentumBasis>,
        subspace: &'a SubspaceData,
        probe_rotation: &'a Rotation,
    ) -> Self {
        let e_cm = F::from_f64(settings.kinematics.e_cm);

        Self {
            real_mass_vector: masses,
            e_cm,
            graph,
            settings,
            esurface_collection,
            overlap_structure,
            transformed_kinematic_point,
            all_lmbs,
            subspace,
            probe_rotation,
        }
    }

    fn new_overlap_builder(&'a self, overlap_group: &'a OverlapGroup) -> OverlapBuilder<'a, T> {
        let subspace = self.subspace;

        let center = overlap_group.center.cast();

        let shifted_loop_momenta = self
            .transformed_kinematic_point
            .representative_sample()
            .loop_moms()
            - &center;

        let radius = shifted_loop_momenta
            .hyper_radius_squared(Some(&subspace.iter_lmb_indices().collect_vec()))
            .sqrt();
        let unit_shifted_momenta = shifted_loop_momenta.rescale(
            &radius.inv(),
            Some(&subspace.iter_lmb_indices().collect_vec()),
        );

        OverlapBuilder {
            counterterm_builder: self,
            overlap_group,
            center,
            unit_shifted_momenta,
            radius,
        }
    }
}

struct OverlapBuilder<'a, T: FloatLike> {
    counterterm_builder: &'a CounterTermBuilder<'a, T>,
    overlap_group: &'a OverlapGroup,
    /// Solver-derived centers already belong to the current probe and cut-side LMB frame.
    center: LoopMomenta<F<T>>,
    unit_shifted_momenta: LoopMomenta<F<T>>,
    radius: F<T>,
}

impl<'a, T: FloatLike> OverlapBuilder<'a, T> {
    fn new_esurface_builder(
        &'a self,
        existing_esurface_id: ExistingEsurfaceId,
    ) -> EsurfaceCTBuilder<'a, T> {
        let esurface_id = self
            .counterterm_builder
            .overlap_structure
            .existing_esurfaces[existing_esurface_id];

        EsurfaceCTBuilder {
            overlap_builder: self,
            _existing_esurface_id: existing_esurface_id,
            esurface: &self.counterterm_builder.esurface_collection[esurface_id],
            esurface_id,
        }
    }
}

const MAX_ITERATIONS: usize = 40;

struct EsurfaceCTBuilder<'a, T: FloatLike> {
    overlap_builder: &'a OverlapBuilder<'a, T>,
    _existing_esurface_id: ExistingEsurfaceId,
    esurface: &'a Esurface,
    esurface_id: EsurfaceID,
}

impl<'a, T: FloatLike> EsurfaceCTBuilder<'a, T> {
    fn solve_rstar(
        self,
        rstar_t_dependence_evaluator: &mut RstarTDependenceEvaluator,
        radial_root_identity: &RadialRootIdentity,
        radial_root_diagnostics: &mut RadialRootDiagnostics,
    ) -> Option<RstarSolution<'a, T>> {
        let subspace = self.overlap_builder.counterterm_builder.subspace;
        let graph = self.overlap_builder.counterterm_builder.graph;
        let lmbs = self.overlap_builder.counterterm_builder.all_lmbs;
        let masses = self.overlap_builder.counterterm_builder.real_mass_vector;

        debug!("subspace: {:?}", subspace);

        let representative_sample = self
            .overlap_builder
            .counterterm_builder
            .transformed_kinematic_point
            .representative_sample();
        let mut center_with_fixed_complement = representative_sample.loop_moms().clone();
        for loop_index in subspace.iter_lmb_indices() {
            center_with_fixed_complement[loop_index] =
                self.overlap_builder.center[loop_index].clone();
        }

        let center_surface_values = self
            .overlap_builder
            .overlap_group
            .existing_esurfaces
            .iter()
            .map(|&existing_esurface_id| {
                let esurface_id = self
                    .overlap_builder
                    .counterterm_builder
                    .overlap_structure
                    .existing_esurfaces[existing_esurface_id];
                let esurface =
                    &self.overlap_builder.counterterm_builder.esurface_collection[esurface_id];
                let value = esurface.compute_from_momenta(
                    subspace.get_lmb(lmbs),
                    masses,
                    &center_with_fixed_complement,
                    representative_sample.external_moms(),
                );
                let is_valid = esurface_value_is_strictly_inside(
                    &value,
                    &self.overlap_builder.counterterm_builder.e_cm,
                );
                (esurface_id, value, is_valid)
            })
            .collect_vec();
        let all_center_values_valid = center_surface_values
            .iter()
            .all(|(_, _, is_valid)| *is_valid);

        crate::debug_tags!(#integration, #subtraction, #threshold, #inspect, #center;
            stage = "lu_threshold_center_values",
            graph = %graph.name,
            selected_esurface_id = self.esurface_id.0,
            rotation_id = %self.overlap_builder.counterterm_builder.probe_rotation.method,
            center_provenance = "current_probe_cut_lmb_frame",
            subspace_loop_indices = ?subspace.iter_lmb_indices().collect_vec(),
            file.active_center = %format!("{}", self.overlap_builder.center),
            file.center_with_fixed_complement = %format!("{}", center_with_fixed_complement),
            file.surface_values = %center_surface_values
                .iter()
                .map(|(esurface_id, value, is_valid)| format!(
                    "local={} value={:+16e} inside={}",
                    esurface_id.0,
                    value,
                    is_valid,
                ))
                .join("\n"),
            "LU threshold center values"
        );

        if !all_center_values_valid {
            warn!(
                graph = %graph.name,
                selected_esurface_id = self.esurface_id.0,
                rotation_id = %self.overlap_builder.counterterm_builder.probe_rotation.method,
                center_provenance = "current_probe_cut_lmb_frame",
                center = %center_with_fixed_complement,
                surface_values = %center_surface_values
                    .iter()
                    .map(|(esurface_id, value, is_valid)| format!(
                        "local={} value={:+16e} inside={}",
                        esurface_id.0,
                        value,
                        is_valid,
                    ))
                    .join("; "),
                "refusing to evaluate an LU threshold counterterm with an invalid probe-frame overlap center"
            );
            return None;
        }

        let (raw_radius_guess, _) = self.esurface.get_radius_guess_subspace(
            &self.overlap_builder.unit_shifted_momenta,
            self.overlap_builder
                .counterterm_builder
                .transformed_kinematic_point
                .representative_sample()
                .external_moms(),
            subspace,
            lmbs,
            graph,
            masses,
        );

        let function = |r: &_| {
            self.esurface.compute_self_and_r_derivative_subspace(
                r,
                &self.overlap_builder.unit_shifted_momenta,
                &self.overlap_builder.center,
                self.overlap_builder
                    .counterterm_builder
                    .transformed_kinematic_point
                    .representative_sample()
                    .external_moms(),
                self.overlap_builder.counterterm_builder.real_mass_vector,
                subspace,
                lmbs,
                graph,
            )
        };

        let zero = raw_radius_guess.zero();
        let mut radius_guess = raw_radius_guess.clone();
        if radius_guess.is_nan() || radius_guess.is_infinite() || radius_guess <= zero {
            radius_guess = self.overlap_builder.counterterm_builder.e_cm.clone();
        }

        debug!("initial radius guess: {:?}", radius_guess);

        // Some residual is expected when Newton stagnates at the representable value nearest the
        // root. Direct convergence uses the active precision; the shared diagnostics may also
        // certify a residual-limited higher-precision result against the preceding precision.
        let tolerance = F::from_f64(
            self.overlap_builder
                .counterterm_builder
                .settings
                .subtraction
                .radial_root_residual_tolerance,
        );
        let solution = match radial_root_diagnostics.solve(
            radial_root_identity,
            &zero,
            &radius_guess,
            function,
            &tolerance,
            MAX_ITERATIONS,
            64,
            &self.overlap_builder.counterterm_builder.e_cm,
        ) {
            Ok(solution) => solution,
            Err(error) => {
                warn!(
                    graph = %graph.name,
                    esurface_id = self.esurface_id.0,
                    rotation_id = %self.overlap_builder.counterterm_builder.probe_rotation.method,
                    center_provenance = "current_probe_cut_lmb_frame",
                    raw_radius_guess = %raw_radius_guess,
                    radius_guess = %radius_guess,
                    error = %error,
                    "refusing to evaluate a threshold counterterm with an invalid radial solution"
                );
                return None;
            }
        };

        debug!("r* solution: {:?}", solution);

        let t_dependent_solution = if rstar_t_dependence_evaluator.supports_t_derivatives() {
            let t_star = match &self
                .overlap_builder
                .counterterm_builder
                .transformed_kinematic_point
                .non_dual_cut_params()
                .tstar
            {
                DualOrNot::NonDual(t_star) => t_star.clone(),
                DualOrNot::Dual(_) => {
                    unreachable!("representative LU cut parameters must stay non-dual")
                }
            };

            Some(
                rstar_t_dependence_evaluator.evaluate(RstarTDependenceInput {
                    t_star: &t_star,
                    radius_star: &solution.solution,
                    overlap_center: &self.overlap_builder.center,
                    subspace,
                    unrescaled_momentum_sample: self
                        .overlap_builder
                        .counterterm_builder
                        .transformed_kinematic_point
                        .unrescaled_sample(),
                    masses,
                    threshold_esurface: self.esurface,
                    lmb: &self
                        .overlap_builder
                        .counterterm_builder
                        .graph
                        .loop_momentum_basis,
                    all_lmbs: lmbs,
                }),
            )
        } else {
            None
        };

        Some(RstarSolution {
            esurface_ct_builder: self,
            solution,
            t_dependent_solution,
        })
    }
}

struct RstarSolution<'a, T: FloatLike> {
    esurface_ct_builder: EsurfaceCTBuilder<'a, T>,
    solution: NewtonIterationResult<T>,
    t_dependent_solution: Option<HyperDual<F<T>>>,
}

struct DualRstarGeometry<T: FloatLike> {
    radius: HyperDual<F<T>>,
    radius_star: HyperDual<F<T>>,
    esurface_derivative: HyperDual<F<T>>,
    rstar_loop_momenta: LoopMomenta<HyperDual<F<T>>>,
    external_moms: ExternalFourMomenta<HyperDual<F<T>>>,
}

impl<'a, T: FloatLike> RstarSolution<'a, T> {
    fn max_supported_order(&self) -> usize {
        self.t_dependent_solution
            .as_ref()
            .map(|dual| dual.values.len().saturating_sub(1))
            .unwrap_or(0)
    }

    fn base_rstar_loop_momenta(&self) -> LoopMomenta<F<T>> {
        let subspace = self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .subspace;

        &self
            .esurface_ct_builder
            .overlap_builder
            .unit_shifted_momenta
            .rescale(&self.solution.solution, subspace.as_subspace_simple())
            + &self.esurface_ct_builder.overlap_builder.center
    }

    fn truncated_rstar_solution(&self, order: usize) -> HyperDual<F<T>> {
        let dual_solution = self
            .t_dependent_solution
            .as_ref()
            .expect("higher-order LU threshold evaluation requires cached r_star(t)");
        HyperDual::from_values(
            simple_n_deriv_shape(order),
            dual_solution.values[..=order].to_vec(),
        )
    }

    fn embedded_truncated_rstar_solution(&self, cut_cff_index: &CutCFFIndex) -> HyperDual<F<T>> {
        let lu_order = cut_cff_index
            .lu_cut_order
            .expect("mixed LU geometry requires lu_cut_order")
            - 1;
        let target_shape = shape_from_cut_cff_index(cut_cff_index)
            .map(HyperDual::new)
            .expect("mixed LU geometry requires a dual shape");
        let t_variable = variable_indices_from_cut_cff_index(cut_cff_index).lu_cut;
        let t_dual = if lu_order == 0 {
            HyperDual::from_values(
                simple_n_deriv_shape(0),
                vec![self.solution.solution.clone()],
            )
        } else {
            self.truncated_rstar_solution(lu_order)
        };

        embed_t_dual_in_target_shape(&t_dual, &target_shape, t_variable)
    }

    fn dual_geometry_for_order(&self, order: usize) -> DualRstarGeometry<T> {
        let source_sample = self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .transformed_kinematic_point
            .sample_for_order(order);
        let dual_loop_momenta = source_sample
            .sample
            .dual_loop_moms
            .as_ref()
            .expect("higher-order LU threshold evaluation requires dual loop momenta")
            .clone();

        let radius_star = self.truncated_rstar_solution(order);
        let center = dualize_loop_momenta(
            &radius_star,
            &self.esurface_ct_builder.overlap_builder.center,
        );
        let external_moms = dualize_external_momenta(&radius_star, source_sample.external_moms());

        let shifted_momenta = dual_loop_momenta
            .iter()
            .zip(center.iter())
            .map(|(momentum, center)| momentum.clone() - center.clone())
            .collect::<LoopMomenta<_>>();

        let radius = dual_shifted_radius(
            &shifted_momenta,
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .subspace,
        );
        let inverse_radius = new_constant(&radius, &radius.values[0].one()) / radius.clone();
        let unit_shifted_momenta = shifted_momenta.rescale(
            &inverse_radius,
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .subspace
                .as_subspace_simple(),
        );

        let rstar_loop_momenta = unit_shifted_momenta
            .rescale(
                &radius_star,
                self.esurface_ct_builder
                    .overlap_builder
                    .counterterm_builder
                    .subspace
                    .as_subspace_simple(),
            )
            .iter()
            .zip(center.iter())
            .map(|(momentum, center)| momentum.clone() + center.clone())
            .collect();

        let (_, esurface_derivative) = compute_self_and_r_derivative_subspace_dual(
            self.esurface_ct_builder.esurface,
            &radius_star,
            &unit_shifted_momenta,
            &center,
            &external_moms,
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .real_mass_vector,
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .subspace,
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .all_lmbs,
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .graph,
        );

        DualRstarGeometry {
            radius,
            radius_star,
            esurface_derivative,
            rstar_loop_momenta,
            external_moms,
        }
    }

    fn dual_geometry_for_cut_cff_index(
        &self,
        cut_cff_index: &CutCFFIndex,
        threshold_variable: Option<usize>,
    ) -> DualRstarGeometry<T> {
        let lu_order = cut_cff_index
            .lu_cut_order
            .expect("mixed LU geometry requires lu_cut_order")
            - 1;
        let source_sample = self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .transformed_kinematic_point
            .sample_for_order(lu_order);
        let target_shape = shape_from_cut_cff_index(cut_cff_index)
            .map(HyperDual::new)
            .expect("mixed LU geometry requires a dual shape");
        let t_variable = variable_indices_from_cut_cff_index(cut_cff_index).lu_cut;
        let dual_loop_momenta = embedded_dual_loop_momenta_for_cut_cff_index(
            source_sample,
            &Some(target_shape.clone()),
            t_variable,
        )
        .expect("mixed LU geometry requires dual loop momenta");

        let mut radius_star = self.embedded_truncated_rstar_solution(cut_cff_index);
        activate_threshold_variable_in_target_shape(&mut radius_star, threshold_variable);
        let center = dualize_loop_momenta(
            &radius_star,
            &self.esurface_ct_builder.overlap_builder.center,
        );
        let external_moms = dualize_external_momenta(&radius_star, source_sample.external_moms());

        let shifted_momenta = dual_loop_momenta
            .iter()
            .zip(center.iter())
            .map(|(momentum, center)| momentum.clone() - center.clone())
            .collect::<LoopMomenta<_>>();

        let radius = dual_shifted_radius(
            &shifted_momenta,
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .subspace,
        );
        let inverse_radius = new_constant(&radius, &radius.values[0].one()) / radius.clone();
        let unit_shifted_momenta = shifted_momenta.rescale(
            &inverse_radius,
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .subspace
                .as_subspace_simple(),
        );

        let rstar_loop_momenta = unit_shifted_momenta
            .rescale(
                &radius_star,
                self.esurface_ct_builder
                    .overlap_builder
                    .counterterm_builder
                    .subspace
                    .as_subspace_simple(),
            )
            .iter()
            .zip(center.iter())
            .map(|(momentum, center)| momentum.clone() + center.clone())
            .collect();

        let (_, esurface_derivative) = compute_self_and_r_derivative_subspace_dual(
            self.esurface_ct_builder.esurface,
            &radius_star,
            &unit_shifted_momenta,
            &center,
            &external_moms,
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .real_mass_vector,
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .subspace,
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .all_lmbs,
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .graph,
        );

        DualRstarGeometry {
            radius,
            radius_star,
            esurface_derivative,
            rstar_loop_momenta,
            external_moms,
        }
    }

    fn non_dual_multichanneling_factor(&self, rstar_sample: &MomentumSample<T>) -> Complex<F<T>> {
        let subspace = self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .subspace;

        let multi_channeling_denominator = self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .overlap_structure
            .overlap_groups
            .iter()
            .map(|group| {
                group
                    .complement
                    .iter()
                    .map(|existing_esurface_id| {
                        let esurface_id = self
                            .esurface_ct_builder
                            .overlap_builder
                            .counterterm_builder
                            .overlap_structure
                            .existing_esurfaces[*existing_esurface_id];
                        let esurface = &self
                            .esurface_ct_builder
                            .overlap_builder
                            .counterterm_builder
                            .esurface_collection[esurface_id];
                        let esurface_value = esurface.compute_from_momenta(
                            subspace.get_lmb(
                                self.esurface_ct_builder
                                    .overlap_builder
                                    .counterterm_builder
                                    .all_lmbs,
                            ),
                            self.esurface_ct_builder
                                .overlap_builder
                                .counterterm_builder
                                .real_mass_vector,
                            rstar_sample.loop_moms(),
                            rstar_sample.external_moms(),
                        );

                        &esurface_value * &esurface_value
                    })
                    .fold(rstar_sample.one(), |acc, value| acc * value)
            })
            .fold(rstar_sample.zero(), |acc, value| acc + value);

        let multichanneling_numerator = self
            .esurface_ct_builder
            .overlap_builder
            .overlap_group
            .complement
            .iter()
            .map(|existing_esurface_id| {
                let esurface_id = self
                    .esurface_ct_builder
                    .overlap_builder
                    .counterterm_builder
                    .overlap_structure
                    .existing_esurfaces[*existing_esurface_id];
                let esurface = &self
                    .esurface_ct_builder
                    .overlap_builder
                    .counterterm_builder
                    .esurface_collection[esurface_id];
                let esurface_value = esurface.compute_from_momenta(
                    subspace.get_lmb(
                        self.esurface_ct_builder
                            .overlap_builder
                            .counterterm_builder
                            .all_lmbs,
                    ),
                    self.esurface_ct_builder
                        .overlap_builder
                        .counterterm_builder
                        .real_mass_vector,
                    rstar_sample.loop_moms(),
                    rstar_sample.external_moms(),
                );

                &esurface_value * &esurface_value
            })
            .fold(rstar_sample.one(), |acc, value| acc * value);

        Complex::new_re(multichanneling_numerator / multi_channeling_denominator)
    }

    fn dual_multichanneling_factor(&self, geometry: &DualRstarGeometry<T>) -> HyperDual<F<T>> {
        let subspace = self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .subspace;
        let lmb = subspace.get_lmb(
            self.esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .all_lmbs,
        );
        let zero = new_constant(&geometry.radius, &geometry.radius.values[0].zero());
        let one = new_constant(&geometry.radius, &geometry.radius.values[0].one());

        let multi_channeling_denominator = self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .overlap_structure
            .overlap_groups
            .iter()
            .map(|group| {
                group
                    .complement
                    .iter()
                    .map(|existing_esurface_id| {
                        let esurface_id = self
                            .esurface_ct_builder
                            .overlap_builder
                            .counterterm_builder
                            .overlap_structure
                            .existing_esurfaces[*existing_esurface_id];
                        let esurface = &self
                            .esurface_ct_builder
                            .overlap_builder
                            .counterterm_builder
                            .esurface_collection[esurface_id];
                        let esurface_value = esurface.compute_from_dual_momenta(
                            lmb,
                            self.esurface_ct_builder
                                .overlap_builder
                                .counterterm_builder
                                .real_mass_vector,
                            &geometry.rstar_loop_momenta,
                            &geometry.external_moms,
                        );

                        esurface_value.clone() * esurface_value
                    })
                    .fold(one.clone(), |acc, value| acc * value)
            })
            .fold(zero.clone(), |acc, value| acc + value);

        let multi_channeling_numerator = self
            .esurface_ct_builder
            .overlap_builder
            .overlap_group
            .complement
            .iter()
            .map(|existing_esurface_id| {
                let esurface_id = self
                    .esurface_ct_builder
                    .overlap_builder
                    .counterterm_builder
                    .overlap_structure
                    .existing_esurfaces[*existing_esurface_id];
                let esurface = &self
                    .esurface_ct_builder
                    .overlap_builder
                    .counterterm_builder
                    .esurface_collection[esurface_id];
                let esurface_value = esurface.compute_from_dual_momenta(
                    lmb,
                    self.esurface_ct_builder
                        .overlap_builder
                        .counterterm_builder
                        .real_mass_vector,
                    &geometry.rstar_loop_momenta,
                    &geometry.external_moms,
                );

                esurface_value.clone() * esurface_value
            })
            .fold(one, |acc, value| acc * value);

        multi_channeling_numerator / multi_channeling_denominator
    }

    fn rstar_sample_for_order<'solution>(
        &'solution self,
        order: usize,
    ) -> RstarSample<'solution, 'a, T> {
        let base_rstar_loop_momenta = self.base_rstar_loop_momenta();

        if order == 0 {
            let mut rstar_sample = self
                .esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .transformed_kinematic_point
                .representative_sample()
                .clone();
            rstar_sample.sample.loop_moms = base_rstar_loop_momenta;
            rstar_sample.sample.dual_loop_moms = None;

            let radius = self.esurface_ct_builder.overlap_builder.radius.clone();
            let radius_star = self.solution.solution.clone();
            let esurface_derivative = self.solution.derivative_at_solution.clone();
            let e_cm = &self
                .esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .e_cm;
            let uv_localisation_settings = &self
                .esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .settings
                .subtraction
                .local_ct_settings
                .uv_localisation;

            let threshold_params = ThresholdParams {
                radius: DualOrNot::NonDual(radius.clone()),
                radius_star: DualOrNot::NonDual(radius_star.clone()),
                esurface_derivative: DualOrNot::NonDual(esurface_derivative.clone()),
                uv_damp_plus: DualOrNot::NonDual(evaluate_uv_damper(
                    &radius,
                    &radius_star,
                    e_cm,
                    uv_localisation_settings,
                )),
                uv_damp_minus: DualOrNot::NonDual(evaluate_uv_damper(
                    &-radius.clone(),
                    &radius_star,
                    e_cm,
                    uv_localisation_settings,
                )),
                h_function: DualOrNot::NonDual(evaluate_integrated_ct_normalisation(
                    &radius,
                    &radius_star,
                    e_cm,
                    &self
                        .esurface_ct_builder
                        .overlap_builder
                        .counterterm_builder
                        .settings
                        .subtraction
                        .integrated_ct_settings,
                )),
            };

            let value_of_multi_channeling_factor =
                DualOrNot::NonDual(self.non_dual_multichanneling_factor(&rstar_sample));

            return RstarSample {
                rstar_solution: self,
                rstar_sample,
                threshold_params,
                value_of_multi_channeling_factor,
            };
        }

        assert!(
            order <= self.max_supported_order(),
            "requested LU threshold derivative order {} but only {} orders are cached",
            order,
            self.max_supported_order()
        );

        let geometry = self.dual_geometry_for_order(order);
        let mut rstar_sample = self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .transformed_kinematic_point
            .sample_for_order(order)
            .clone();
        rstar_sample.sample.loop_moms = base_rstar_loop_momenta;
        rstar_sample.sample.dual_loop_moms = Some(geometry.rstar_loop_momenta.clone());

        let e_cm = &self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .e_cm;
        let uv_localisation_settings = &self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .settings
            .subtraction
            .local_ct_settings
            .uv_localisation;

        let threshold_params = ThresholdParams {
            radius: DualOrNot::Dual(geometry.radius.clone()),
            radius_star: DualOrNot::Dual(geometry.radius_star.clone()),
            esurface_derivative: DualOrNot::Dual(geometry.esurface_derivative.clone()),
            uv_damp_plus: DualOrNot::Dual(evaluate_uv_damper_dual(
                &geometry.radius,
                &geometry.radius_star,
                e_cm,
                uv_localisation_settings,
            )),
            uv_damp_minus: DualOrNot::Dual(evaluate_uv_damper_dual(
                &-geometry.radius.clone(),
                &geometry.radius_star,
                e_cm,
                uv_localisation_settings,
            )),
            h_function: DualOrNot::Dual(evaluate_integrated_ct_normalisation_dual(
                &geometry.radius,
                &geometry.radius_star,
                e_cm,
                &self
                    .esurface_ct_builder
                    .overlap_builder
                    .counterterm_builder
                    .settings
                    .subtraction
                    .integrated_ct_settings,
            )),
        };

        let value_of_multi_channeling_factor = DualOrNot::Dual(real_dual_to_complex(
            self.dual_multichanneling_factor(&geometry),
        ));

        RstarSample {
            rstar_solution: self,
            rstar_sample,
            threshold_params,
            value_of_multi_channeling_factor,
        }
    }

    fn rstar_sample_for_cut_cff_index<'solution>(
        &'solution self,
        cut_cff_index: &CutCFFIndex,
        threshold_variable: Option<usize>,
    ) -> RstarSample<'solution, 'a, T> {
        let lu_order = cut_cff_index
            .lu_cut_order
            .expect("LU threshold sampling requires lu_cut_order")
            - 1;

        if shape_from_cut_cff_index(cut_cff_index).is_none() {
            return self.rstar_sample_for_order(lu_order);
        }

        let base_rstar_loop_momenta = self.base_rstar_loop_momenta();
        let geometry = self.dual_geometry_for_cut_cff_index(cut_cff_index, threshold_variable);
        let mut rstar_sample = self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .transformed_kinematic_point
            .sample_for_order(lu_order)
            .clone();
        rstar_sample.sample.loop_moms = base_rstar_loop_momenta;
        rstar_sample.sample.dual_loop_moms = Some(geometry.rstar_loop_momenta.clone());

        let e_cm = &self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .e_cm;
        let uv_localisation_settings = &self
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .settings
            .subtraction
            .local_ct_settings
            .uv_localisation;

        let threshold_params = ThresholdParams {
            radius: DualOrNot::Dual(geometry.radius.clone()),
            radius_star: DualOrNot::Dual(geometry.radius_star.clone()),
            esurface_derivative: DualOrNot::Dual(geometry.esurface_derivative.clone()),
            uv_damp_plus: DualOrNot::Dual(evaluate_uv_damper_dual(
                &geometry.radius,
                &geometry.radius_star,
                e_cm,
                uv_localisation_settings,
            )),
            uv_damp_minus: DualOrNot::Dual(evaluate_uv_damper_dual(
                &-geometry.radius.clone(),
                &geometry.radius_star,
                e_cm,
                uv_localisation_settings,
            )),
            h_function: DualOrNot::Dual(evaluate_integrated_ct_normalisation_dual(
                &geometry.radius,
                &geometry.radius_star,
                e_cm,
                &self
                    .esurface_ct_builder
                    .overlap_builder
                    .counterterm_builder
                    .settings
                    .subtraction
                    .integrated_ct_settings,
            )),
        };

        let value_of_multi_channeling_factor = DualOrNot::Dual(real_dual_to_complex(
            self.dual_multichanneling_factor(&geometry),
        ));

        RstarSample {
            rstar_solution: self,
            rstar_sample,
            threshold_params,
            value_of_multi_channeling_factor,
        }
    }
}

struct RstarSample<'solution, 'a, T: FloatLike> {
    rstar_solution: &'solution RstarSolution<'a, T>,
    rstar_sample: MomentumSample<T>,
    threshold_params: ThresholdParams<T>,
    value_of_multi_channeling_factor: DualOrNot<Complex<F<T>>>,
}

impl<'solution, 'a, T: FloatLike> RstarSample<'solution, 'a, T> {
    fn extract_threshold_parameters(&self, is_first_call: bool) -> ThresholdParams<T> {
        if is_first_call {
            let edges_in_esurface = self
                .rstar_solution
                .esurface_ct_builder
                .esurface
                .energies
                .iter()
                .map(|i| {
                    self.rstar_solution
                        .esurface_ct_builder
                        .overlap_builder
                        .counterterm_builder
                        .graph[i]
                        .0
                        .name
                        .clone()
                })
                .collect_vec();

            debug!("esurface_id: {}", self.get_esurface_id().0);
            debug!("edges in esurface: {:?}", edges_in_esurface);
            debug!("radius: {}", self.threshold_params.radius);
            debug!("radius_star: {}", self.threshold_params.radius_star);
            debug!(
                "esurface_derivative: {}",
                self.threshold_params.esurface_derivative
            );
            debug!("uv_damp_plus: {}", self.threshold_params.uv_damp_plus);
            debug!("uv_damp_minus: {}", self.threshold_params.uv_damp_minus);
            debug!("h_function: {}", self.threshold_params.h_function);
            debug!(
                "value of multi-channeling factor: {}",
                self.value_of_multi_channeling_factor
            );
        }

        self.threshold_params.clone()
    }

    fn get_inverse_transformed_sample(&self) -> MomentumSample<T> {
        let subspace = self
            .rstar_solution
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .subspace;
        let current_lmb = subspace.get_lmb(
            self.rstar_solution
                .esurface_ct_builder
                .overlap_builder
                .counterterm_builder
                .all_lmbs,
        );
        let target_lmb = &self
            .rstar_solution
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .graph
            .loop_momentum_basis;

        self.rstar_sample.lmb_transform(current_lmb, target_lmb)
    }

    fn get_esurface_id(&self) -> EsurfaceID {
        self.rstar_solution.esurface_ct_builder.esurface_id
    }
}

fn merge_and_inverse_transform<T: FloatLike>(
    left_sample: &RstarSample<'_, '_, T>,
    right_sample: &RstarSample<'_, '_, T>,
) -> MomentumSample<T> {
    let left_subspace = left_sample
        .rstar_solution
        .esurface_ct_builder
        .overlap_builder
        .counterterm_builder
        .subspace;

    let right_subspace = right_sample
        .rstar_solution
        .esurface_ct_builder
        .overlap_builder
        .counterterm_builder
        .subspace;

    assert!(
        left_subspace.is_mergable_with(right_subspace),
        "incompatible subspaces for merging samples"
    );

    let mut merged_sample = left_sample.rstar_sample.clone();
    for lmb_index in right_subspace.iter_lmb_indices() {
        let right_momentum = right_sample.rstar_sample.loop_moms()[lmb_index].clone();
        merged_sample.sample.loop_moms[lmb_index] = right_momentum;
    }

    match (
        merged_sample.sample.dual_loop_moms.as_mut(),
        right_sample.rstar_sample.sample.dual_loop_moms.as_ref(),
    ) {
        (Some(merged_dual_loop_moms), Some(right_dual_loop_moms)) => {
            for lmb_index in right_subspace.iter_lmb_indices() {
                merged_dual_loop_moms[lmb_index] = right_dual_loop_moms[lmb_index].clone();
            }
        }
        (None, None) => {}
        _ => {
            unreachable!("iterated LU samples must either both carry dual loop momenta or neither")
        }
    }

    let current_lmb = left_subspace.get_lmb(
        left_sample
            .rstar_solution
            .esurface_ct_builder
            .overlap_builder
            .counterterm_builder
            .all_lmbs,
    );

    let target_lmb = &left_sample
        .rstar_solution
        .esurface_ct_builder
        .overlap_builder
        .counterterm_builder
        .graph
        .loop_momentum_basis;

    merged_sample.lmb_transform(current_lmb, target_lmb)
}

#[cfg(test)]
mod tests {
    use super::{LUCounterTerm, extract_zero_threshold_coefficient_t_dual};
    use crate::cff::CutCFFIndex;
    use crate::processes::CutGroupId;
    use crate::utils::{
        F,
        hyperdual_utils::{new_from_values, shape_from_cut_cff_index},
    };
    use symbolica::domains::dual::HyperDual;
    use typed_index_collections::{TiVec, ti_vec};

    #[test]
    fn inactive_cut_group_guard_reports_runtime_access() {
        let counterterm = LUCounterTerm {
            evaluators: TiVec::new(),
            thresholds: TiVec::new(),
            subspaces: TiVec::new(),
            rstar_dependence_calculator: TiVec::new(),
            active_cut_groups: ti_vec![false],
            active_left_thresholds: ti_vec![TiVec::new()],
            active_right_thresholds: ti_vec![TiVec::new()],
            active_iterated_thresholds: TiVec::new(),
        };

        let error = counterterm
            .ensure_active_cut_group(CutGroupId(0))
            .unwrap_err();
        assert!(error.to_string().contains("generation marked it inactive"));
    }

    #[test]
    fn zero_threshold_coefficient_lookup_is_explicit() {
        let shape = HyperDual::new(
            shape_from_cut_cff_index(&CutCFFIndex {
                left_threshold_order: Some(2),
                right_threshold_order: Some(2),
                lu_cut_order: Some(2),
            })
            .unwrap(),
        );
        let mixed_dual = new_from_values(
            &shape,
            &[
                F(10.0_f64),
                F(11.0_f64),
                F(20.0_f64),
                F(30.0_f64),
                F(40.0_f64),
                F(21.0_f64),
                F(31.0_f64),
                F(41.0_f64),
            ],
        );

        let zero_threshold = extract_zero_threshold_coefficient_t_dual(&mixed_dual, 0);

        assert_eq!(zero_threshold.values, vec![F(10.0_f64), F(11.0_f64)]);
    }
}