gammalooprs/cff/
expression.rs

1use std::collections::{BTreeMap, BTreeSet};
2
3pub use three_dimensional_reps::expression::{
4    AllOrientations, CFFVariant, OrientationData, OrientationExpression, OrientationID,
5    OrientationSelector, RaisedEsurfaceData, RaisedEsurfaceDataView, RaisedEsurfaceGroup,
6    RaisedEsurfaceGroupView, RaisedEsurfaceId,
7};
8
9use itertools::Itertools;
10use linnet::half_edge::involution::EdgeIndex;
11use spenso::structure::{
12    abstract_index::AIND_SYMBOLS,
13    representation::{LibraryRep, Minkowski, RepName},
14};
15use symbolica::{
16    atom::{Atom, AtomCore, FunctionBuilder},
17    function,
18    id::Replacement,
19};
20use three_dimensional_reps::MediumMode;
21
22use super::{
23    esurface::Esurface, hsurface::Hsurface, orientations::GraphOrientation,
24    surface::GammaLoopLinearEnergyExpr,
25};
26use crate::{
27    graph::Graph,
28    utils::{GS, W_, ose_atom_from_index},
29};
30
31pub type ThreeDExpression<O> =
32    three_dimensional_reps::expression::ThreeDExpression<O, Esurface, Hsurface>;
33pub type CFFExpression<O> =
34    three_dimensional_reps::expression::CFFExpression<O, Esurface, Hsurface>;
35
36pub(crate) fn normalize_three_d_expression_cut_support_with_raised_edge_groups<O>(
37    expression: &mut ThreeDExpression<O>,
38    raised_edge_groups: &[Vec<EdgeIndex>],
39) where
40    O: From<usize> + Into<usize>,
41{
42    for orientation in expression.orientations.iter_mut() {
43        for variant in &mut orientation.variants {
44            variant.denominator_edges = normalize_cut_edge_support_with_raised_edge_groups(
45                &variant.denominator_edges,
46                raised_edge_groups,
47            );
48            variant.denominator_edge_support_signs = normalize_cut_edge_support_signs(
49                std::mem::take(&mut variant.denominator_edge_support_signs),
50                raised_edge_groups,
51            );
52        }
53    }
54}
55
56fn normalize_cut_edge_support_signs(
57    support_signs: BTreeMap<Vec<EdgeIndex>, i64>,
58    raised_edge_groups: &[Vec<EdgeIndex>],
59) -> BTreeMap<Vec<EdgeIndex>, i64> {
60    support_signs
61        .into_iter()
62        .fold(BTreeMap::new(), |mut normalized, (support, sign)| {
63            let support =
64                normalize_cut_edge_support_with_raised_edge_groups(&support, raised_edge_groups);
65            *normalized.entry(support).or_insert(1) *= sign;
66            normalized
67        })
68}
69
70pub(crate) fn normalize_cut_edge_support_with_raised_edge_groups(
71    edges: &[EdgeIndex],
72    raised_edge_groups: &[Vec<EdgeIndex>],
73) -> Vec<EdgeIndex> {
74    edges
75        .iter()
76        .map(|edge| {
77            raised_edge_groups
78                .iter()
79                .find(|group| group.contains(edge))
80                .and_then(|group| group.first())
81                .copied()
82                .unwrap_or(*edge)
83        })
84        .sorted()
85        .dedup()
86        .collect()
87}
88
89pub trait GammaLoopCFFVariant {
90    fn to_atom_gs(&self) -> Atom;
91}
92
93impl GammaLoopCFFVariant for CFFVariant {
94    fn to_atom_gs(&self) -> Atom {
95        let half_edge_factor = self
96            .half_edges
97            .iter()
98            .map(|edge_id| Atom::num(1) / (Atom::num(2) * ose_atom_from_index(*edge_id)))
99            .reduce(|acc, factor| acc * factor)
100            .unwrap_or_else(|| Atom::num(1));
101        let scale_factor = if self.uniform_scale_power == 0 {
102            Atom::num(1)
103        } else {
104            Atom::num(1)
105                / Atom::var(GS.numerator_sampling_scale).pow(self.uniform_scale_power as i64)
106        };
107
108        let numerator_surface_factor = self
109            .numerator_surfaces
110            .iter()
111            .map(|surface_id| Atom::from(*surface_id))
112            .reduce(|acc, factor| acc * factor)
113            .unwrap_or_else(|| Atom::num(1));
114
115        Atom::num(self.prefactor.clone())
116            * self.thermal_weight.to_atom()
117            * half_edge_factor
118            * scale_factor
119            * numerator_surface_factor
120            * self.denominator.to_atom_inv()
121    }
122}
123
124pub trait GammaLoopOrientationExpression {
125    fn to_atom_gs(&self) -> Atom;
126    fn energy_replacements_gs(&self, graph: &Graph) -> Vec<Replacement>;
127}
128
129impl GammaLoopOrientationExpression for OrientationExpression {
130    fn to_atom_gs(&self) -> Atom {
131        let expression = self
132            .variants
133            .iter()
134            .map(GammaLoopCFFVariant::to_atom_gs)
135            .reduce(|acc, atom| acc + atom)
136            .unwrap_or_else(Atom::new);
137        if self
138            .variants
139            .iter()
140            .any(|variant| variant.thermal_weight.medium_mode != MediumMode::Vacuum)
141        {
142            // Each coefficient must carry its own chemical-potential orientation
143            // before an explicit sum removes the residue-map selectors.
144            self.data.orientation.select(expression)
145        } else {
146            expression
147        }
148    }
149
150    fn energy_replacements_gs(&self, graph: &Graph) -> Vec<Replacement> {
151        energy_map_replacements_gs(
152            self.edge_energy_map
153                .iter()
154                .map(|energy| energy.to_atom_gs(&[])),
155            graph,
156        )
157    }
158}
159
160pub(crate) fn energy_map_replacements_gs(
161    edge_energy_map: impl IntoIterator<Item = Atom>,
162    graph: &Graph,
163) -> Vec<Replacement> {
164    let edge_energy_map = edge_energy_map.into_iter().collect::<Vec<_>>();
165    let mut replacements = Vec::new();
166    let mink_index = LibraryRep::from(Minkowski {}).to_symbolic([Atom::var(W_.a__)]);
167    let external_edges = graph
168        .underlying
169        .iter_edges()
170        .filter_map(|(pair, edge_id, _)| (!pair.is_paired()).then_some(edge_id))
171        .collect::<BTreeSet<_>>();
172
173    for (edge_id, energy) in edge_energy_map.iter().enumerate() {
174        let edge_id = EdgeIndex(edge_id);
175        // Remapping a generated internal-edge map into GammaLoop's physical
176        // namespace pads unpaired external-edge slots with zero. Those slots
177        // are not residue samples: their EMR momenta must remain the original
178        // external momenta carried by the factorized numerator.
179        if external_edges.contains(&edge_id) {
180            continue;
181        }
182        replacements.push(Replacement::new(
183            GS.emr_mom(edge_id, AIND_SYMBOLS.cind.call(Atom::Zero))
184                .to_pattern(),
185            energy.clone().to_pattern(),
186        ));
187        replacements.push(Replacement::new(
188            GS.emr_mom(edge_id, &mink_index).to_pattern(),
189            (GS.emr_vec(edge_id, &mink_index) + energy * GS.energy_delta(&mink_index)).to_pattern(),
190        ));
191    }
192
193    for (loop_id, loop_edge_id) in graph.loop_momentum_basis.loop_edges.iter_enumerated() {
194        let loop_id = usize::from(loop_id);
195        let loop_id_atom = Atom::num(loop_id as i64);
196        let energy = edge_energy_map
197            .get(usize::from(*loop_edge_id))
198            .cloned()
199            .unwrap_or_else(Atom::new);
200        replacements.push(Replacement::new(
201            function!(
202                GS.loop_mom,
203                loop_id_atom.clone(),
204                AIND_SYMBOLS.cind.call(Atom::Zero)
205            )
206            .to_pattern(),
207            energy.clone().to_pattern(),
208        ));
209        for spatial_index in 1..=3 {
210            replacements.push(Replacement::new(
211                function!(
212                    GS.loop_mom,
213                    loop_id_atom.clone(),
214                    AIND_SYMBOLS.cind.call(spatial_index)
215                )
216                .to_pattern(),
217                GS.emr_mom(*loop_edge_id, AIND_SYMBOLS.cind.call(spatial_index))
218                    .to_pattern(),
219            ));
220        }
221        replacements.push(Replacement::new(
222            FunctionBuilder::new(GS.loop_mom)
223                .add_arg(loop_id as i64)
224                .add_arg(mink_index.as_view())
225                .finish()
226                .to_pattern(),
227            (GS.emr_vec(*loop_edge_id, &mink_index) + energy * GS.energy_delta(&mink_index))
228                .to_pattern(),
229        ));
230    }
231
232    replacements
233}
234
235#[cfg(test)]
236mod tests {
237    use linnet::half_edge::involution::{EdgeVec, Orientation};
238    use three_dimensional_reps::thermal::OCCUPATION_FORM;
239
240    use super::*;
241    use crate::{
242        cff::surface::LinearEnergyExpr, dot, graph::parse::from_dot::IntoGraph,
243        initialisation::test_initialise, utils::external_energy_atom_from_index,
244    };
245
246    #[test]
247    fn vacuum_subtraction_reaches_root_and_unit_cographs() -> color_eyre::Result<()> {
248        use crate::{
249            graph::cuts::CutSet, integrands::process::param_builder::ParamBuilderGraph,
250            settings::global::OrientationPattern,
251        };
252        test_initialise()?;
253        let mut graph: Graph = dot!(digraph G {
254            edge [particle="scalar_1"]; node [num=1];
255            v1 -> v2; v1 -> v2; v1 -> v2;
256        }, "scalars")?;
257        for medium in [
258            MediumMode::ThermodynamicEquilibrium,
259            MediumMode::ZeroTemperatureEquilibrium,
260        ] {
261            let mut options = graph.denominator_only_cff_3d_expression_options();
262            options.medium_mode = medium;
263            let production =
264                graph.generate_3d_expression_for_integrand(&[], &None, &options, None)?;
265            let cutset = CutSet::empty(graph.n_hedges());
266            let pattern = OrientationPattern::default();
267            let raw = graph
268                .cff_from_production_expression(&production, &cutset, &pattern)?
269                .expression_with_selectors();
270            let actual = graph
271                .cff_from_production_expression(&production, &cutset, &pattern)?
272                .with_vacuum_subtraction(true)
273                .expression_with_selectors();
274            assert_eq!(raw.iter().count(), actual.iter().count());
275            for ((raw_key, raw), (actual_key, actual)) in raw.iter().zip(actual.iter()) {
276                assert_eq!(raw_key, actual_key);
277                let vacuum = graph.make_thermal_distributions_explicit(
278                    raw,
279                    MediumMode::Vacuum,
280                    graph.iter_edges().map(|(_, edge, _)| edge),
281                )?;
282                assert!(raw.is_zero() || raw.contains_symbol(GS.thermal_weight_wrapper));
283                assert!(actual.is_zero() || actual.contains_symbol(GS.thermal_weight_wrapper));
284                // Unwrap only the diagnostic copy when comparing the two algebraic forms.
285                let difference = (actual - (raw - vacuum))
286                    .replace(function!(GS.thermal_weight_wrapper, W_.a_))
287                    .with(W_.a_)
288                    .replace_multiple(&*OCCUPATION_FORM);
289                assert!(difference.expand().is_zero());
290            }
291            let full = graph.full_filter();
292            let unit = graph
293                .cff(&full, &cutset, &pattern, &options, None)?
294                .with_vacuum_subtraction(true)
295                .expression_with_selectors();
296            assert!(
297                unit.iter().all(|(_, atom)| atom.is_zero()),
298                "the fully contracted cograph must subtract its empty weight"
299            );
300            assert!(
301                !graph
302                    .cff(&full, &cutset, &pattern, &options, None)?
303                    .expression_with_selectors()
304                    .iter()
305                    .all(|(_, atom)| atom.is_zero())
306            );
307        }
308        Ok(())
309    }
310
311    #[test]
312    fn thermal_surface_conversion_preserves_distribution_atoms() -> color_eyre::Result<()> {
313        test_initialise()?;
314        let graph: Graph = dot!(digraph thermal_triangle {
315            edge [num=1 mass=0]
316            node [num=1]
317            ext [style=invis]
318            ext -> a [id=3]
319            ext -> b [id=4]
320            c -> ext [id=5]
321            c -> a [id=0 lmb_id=0]
322            a -> b [id=1]
323            b -> c [id=2]
324        })?;
325        for medium_mode in [
326            MediumMode::ThermodynamicEquilibrium,
327            MediumMode::ZeroTemperatureEquilibrium,
328        ] {
329            let mut graph = graph.clone();
330            let mut options = graph.denominator_only_cff_3d_expression_options();
331            options.medium_mode = medium_mode;
332            let expression = graph
333                .generate_3d_expression_for_integrand(&[], &None, &options, None)?
334                .expression;
335            let atom = expression
336                .orientations
337                .iter()
338                .map(GammaLoopOrientationExpression::to_atom_gs)
339                .fold(Atom::Zero, |sum, orientation| sum + orientation);
340            assert_eq!(atom, expression.to_atom(AllOrientations));
341            assert!(!atom.contains_symbol(GS.sign));
342
343            // Typed weights emit distributions directly; surface conversion must preserve them.
344            let substituted = expression.surfaces.substitute_energies(&atom, &[]);
345            assert!(
346                substituted.contains_symbol(GS.thermal_distribution),
347                "{medium_mode:?} CFF lost its thermal distributions: {substituted}"
348            );
349        }
350        Ok(())
351    }
352
353    #[test]
354    fn affine_energy_map_keeps_two_numerator_factors_under_one_map() -> color_eyre::Result<()> {
355        test_initialise()?;
356        let graph: Graph = dot!(
357            digraph affine_map {
358                edge [num=1 mass=0]
359                node [num=1]
360                A -> B [id=0]
361                A -> B [id=1]
362            }
363        )?;
364        let energy_map = LinearEnergyExpr {
365            internal_terms: vec![(EdgeIndex(0), 2.into()), (EdgeIndex(1), (-3).into())],
366            external_terms: vec![(EdgeIndex(2), 5.into())],
367            uniform_scale_coeff: 7.into(),
368            constant: 11.into(),
369        }
370        .canonical();
371        let orientation = OrientationExpression {
372            data: OrientationData::new(EdgeVec::from_iter([
373                Orientation::Undirected,
374                Orientation::Undirected,
375            ])),
376            loop_energy_map: Vec::new(),
377            edge_energy_map: vec![energy_map.clone(), LinearEnergyExpr::zero()],
378            variants: Vec::new(),
379        };
380        let factor_a = Atom::var(symbolica::symbol!("factor_a"));
381        let factor_b = Atom::var(symbolica::symbol!("factor_b"));
382        let energy_component = GS.emr_mom(EdgeIndex(0), AIND_SYMBOLS.cind.call(Atom::Zero));
383        let numerator =
384            (energy_component.clone() + factor_a.clone()) * (energy_component + factor_b.clone());
385
386        let mapped = numerator.replace_multiple(orientation.energy_replacements_gs(&graph));
387        let mapped_energy = Atom::num(2) * ose_atom_from_index(EdgeIndex(0))
388            - Atom::num(3) * ose_atom_from_index(EdgeIndex(1))
389            + Atom::num(5) * external_energy_atom_from_index(EdgeIndex(2))
390            + Atom::num(7) * Atom::var(GS.numerator_sampling_scale)
391            + Atom::num(11);
392        let expected =
393            (mapped_energy.clone() + factor_a.clone()) * (mapped_energy.clone() + factor_b.clone());
394        assert_eq!(mapped, expected);
395        let mixed = (mapped_energy.clone() + factor_a) * (mapped_energy + Atom::num(1) + factor_b);
396        assert_ne!(mapped, mixed);
397        Ok(())
398    }
399
400    #[test]
401    fn padded_external_energy_slots_do_not_erase_factorized_external_momenta()
402    -> color_eyre::Result<()> {
403        test_initialise()?;
404        let graph: Graph = dot!(digraph padded_external_energy_slots {
405            edge [num=1 mass=0]
406            node [num=1]
407            ext [style=invis]
408            v0;
409            v1;
410            v2;
411            ext -> v0 [id=0]
412            ext -> v1 [id=1]
413            v2 -> ext [id=2]
414            v0 -> v2 [id=3 lmb_id=0]
415            v1 -> v0 [id=4]
416            v2 -> v1 [id=5]
417        })?;
418        let mapped_internal_energy = LinearEnergyExpr {
419            internal_terms: vec![(EdgeIndex(3), 2.into())],
420            external_terms: vec![(EdgeIndex(0), (-1).into())],
421            uniform_scale_coeff: 0.into(),
422            constant: 0.into(),
423        }
424        .canonical();
425        let mut edge_energy_map = vec![LinearEnergyExpr::zero(); graph.underlying.n_edges()];
426        edge_energy_map[3] = mapped_internal_energy.clone();
427        let external_energy = GS.emr_mom(EdgeIndex(0), AIND_SYMBOLS.cind.call(Atom::Zero));
428        let internal_energy = GS.emr_mom(EdgeIndex(3), AIND_SYMBOLS.cind.call(Atom::Zero));
429        let spectator = Atom::var(symbolica::symbol!("spectator"));
430        let numerator = external_energy.clone() * (internal_energy + spectator.clone());
431
432        let mapped = numerator.replace_multiple(energy_map_replacements_gs(
433            edge_energy_map.iter().map(|energy| energy.to_atom_gs(&[])),
434            &graph,
435        ));
436        assert_eq!(
437            mapped,
438            external_energy * (mapped_internal_energy.to_atom_gs(&[]) + spectator)
439        );
440        Ok(())
441    }
442
443    #[test]
444    fn raised_alias_normalization_merges_cut_support_signs() {
445        let support_signs = BTreeMap::from([
446            (vec![EdgeIndex(2), EdgeIndex(3)], -1),
447            (vec![EdgeIndex(3), EdgeIndex(4)], -1),
448        ]);
449
450        assert_eq!(
451            normalize_cut_edge_support_signs(support_signs, &[vec![EdgeIndex(2), EdgeIndex(4)]],),
452            BTreeMap::from([(vec![EdgeIndex(2), EdgeIndex(3)], 1)])
453        );
454    }
455}
gammalooprs/cff/expression.rs

gammalooprs/cff/
expression.rs
