gammalooprs/numerator/
graph.rs

1use linnet::half_edge::{
2    HedgeGraph,
3    involution::{EdgeData, Flow, Orientation},
4    subgraph::{SuBitGraph, SubSetLike},
5};
6use spenso::{
7    iterators::IteratableTensor,
8    network::parsing::ShadowedStructure,
9    shadowing::Concretize,
10    structure::{Canonicalized, TensorStructure},
11    tensors::parametric::ParamTensor,
12};
13use symbolica::atom::{Atom, AtomOrView, FunctionBuilder, Symbol};
14
15use crate::{
16    graph::{Edge, Graph, HedgeData, edge::ParseEdge},
17    model::ArcParticle,
18    utils::GS,
19};
20
21use super::aind::Aind;
22
23pub trait HedgePolarizationData {
24    fn polarization(&self, builder: FunctionBuilder) -> Atom;
25
26    fn edge_spin_slots<'a>(
27        &'a self,
28    ) -> impl Iterator<Item = spenso::structure::representation::LibrarySlot<Aind>> + 'a;
29}
30
31impl HedgePolarizationData for HedgeData {
32    fn polarization(&self, builder: FunctionBuilder) -> Atom {
33        HedgeData::polarization(self, builder)
34    }
35
36    fn edge_spin_slots<'a>(
37        &'a self,
38    ) -> impl Iterator<Item = spenso::structure::representation::LibrarySlot<Aind>> + 'a {
39        HedgeData::edge_spin_slots(self)
40    }
41}
42
43pub trait GeneratePolarizations {
44    /// Returns the polarizations of the given subgraph. One polarization per half-edge.
45    /// If you only want those that are dangling first get the crown of the subgraph.
46    fn generate_polarizations_of<S: SubSetLike>(&self, subgraph: &S) -> Atom;
47    fn generate_polarization_parameters_of<S: SubSetLike>(&self, subgraph: &S) -> Vec<Atom>;
48    // fn generate_polarizations(&self) -> Atom;
49    // fn generate_polarization_params(&self) -> Vec<Atom>;
50}
51
52impl Graph {
53    pub fn generate_polarizations(&self) -> Atom {
54        self.generate_polarizations_of(&self.underlying.external_filter::<SuBitGraph>())
55    }
56    pub fn generate_polarization_params(&self) -> Vec<Atom> {
57        self.generate_polarization_parameters_of(&self.underlying.external_filter::<SuBitGraph>())
58    }
59}
60
61impl GeneratePolarizations for Graph {
62    fn generate_polarizations_of<S: SubSetLike>(&self, subgraph: &S) -> Atom {
63        self.underlying.generate_polarizations_of(subgraph)
64    }
65
66    fn generate_polarization_parameters_of<S: SubSetLike>(&self, subgraph: &S) -> Vec<Atom> {
67        self.underlying
68            .generate_polarization_parameters_of(subgraph)
69    }
70
71    // fn generate_polarizations(&self) -> Atom {
72    //     self.generate_polarizations_of(&self.underlying.external_filter())
73    // }
74    // fn generate_polarization_params(&self) -> Vec<Atom> {
75    //     self.generate_polarization_parameters_of(&self.underlying.external_filter())
76    // }
77
78    // pub(crate) fn polarizations_values(&self) -> Vec<Atom> {
79    //     self.polarizations_of(&self.underlying.external_filter())
80    // }
81}
82
83impl<V, E> GeneratePolarizations for HedgeGraph<E, V, HedgeData>
84where
85    for<'a> EdgeData<&'a E>: ReversibleEdge,
86{
87    fn generate_polarizations_of<S: SubSetLike>(&self, subgraph: &S) -> Atom {
88        let mut pols = Atom::num(1);
89
90        for h in subgraph.included_iter() {
91            let eid = self[&h];
92            pols *= self.get_edge_data_full(h).polarization(
93                &[Atom::num(eid.0)],
94                &self[h],
95                self.flow(h),
96            );
97        }
98
99        pols
100    }
101
102    fn generate_polarization_parameters_of<S: SubSetLike>(&self, subgraph: &S) -> Vec<Atom> {
103        let mut pols = Vec::new();
104
105        for h in subgraph.included_iter() {
106            let eid = self[&h];
107
108            let Some(pol) = self.get_edge_data_full(h).polarization_structure(
109                &[Atom::num(eid.0)],
110                &self[h],
111                self.flow(h),
112            ) else {
113                continue;
114            };
115
116            let (structure, layout) = pol.into_parts();
117            let concrete: ParamTensor<_> = structure
118                .to_shell()
119                .concretize_logical(&layout)
120                .expect("polarization parameters require concrete tensor dimensions");
121
122            for (_, i) in concrete.iter_flat() {
123                pols.push(i.to_owned());
124            }
125        }
126
127        pols
128    }
129
130    // fn generate_polarizations(&self) -> Atom {
131    //     self.generate_polarizations_of(&self.external_filter())
132    // }
133    // fn generate_polarization_params(&self) -> Vec<Atom> {
134    //     self.generate_polarization_parameters_of(&self.external_filter())
135    // }
136
137    // pub(crate) fn polarizations_values(&self) -> Vec<Atom> {
138    //     self.polarizations_of(&self.underlying.external_filter())
139    // }
140}
141
142pub trait ReversibleEdge {
143    fn pdg(&self) -> isize;
144
145    fn pol_symbol(&self, flow: Flow) -> Option<Symbol>;
146
147    fn polarization<'a, T, H>(&self, add_args: &'a [T], hedge_data: &H, flow: Flow) -> Atom
148    where
149        &'a T: Into<AtomOrView<'a>>,
150        H: HedgePolarizationData;
151
152    fn polarization_structure<'a, T, H>(
153        &self,
154        add_args: &'a [T],
155        hedge_data: &H,
156        flow: Flow,
157    ) -> Option<Canonicalized<ShadowedStructure<Aind>>>
158    where
159        &'a T: Into<AtomOrView<'a>>,
160        H: HedgePolarizationData;
161}
162
163impl ReversibleEdge for EdgeData<&Edge> {
164    fn pdg(&self) -> isize {
165        // Parsed edges already store the signed particle species. Their arrow
166        // orientation follows that species and must not conjugate it again.
167        self.data.particle().map_or(0, |particle| particle.pdg_code)
168    }
169
170    fn pol_symbol(&self, flow: Flow) -> Option<Symbol> {
171        if let Some(p) = self.data.particle() {
172            // println!("{}{:?}", p.spin, flow);
173            match (p.spin, flow) {
174                (2, Flow::Sink) => {
175                    if self.pdg() > 0 {
176                        Some(GS.u)
177                    } else {
178                        Some(GS.vbar)
179                    }
180                }
181                (2, Flow::Source) => {
182                    if self.pdg() > 0 {
183                        Some(GS.ubar)
184                    } else {
185                        Some(GS.v)
186                    }
187                }
188                (3, Flow::Sink) => Some(GS.epsilon),
189                (3, Flow::Source) => Some(GS.epsilonbar),
190                _ => None,
191            }
192        } else {
193            None
194        }
195    }
196
197    fn polarization<'a, T, H>(&self, add_args: &'a [T], hedge_data: &H, flow: Flow) -> Atom
198    where
199        &'a T: Into<AtomOrView<'a>>,
200        H: HedgePolarizationData,
201    {
202        if let Some(name) = self.pol_symbol(flow) {
203            hedge_data.polarization(FunctionBuilder::new(name).add_args(add_args))
204        } else {
205            Atom::num(1)
206        }
207    }
208
209    fn polarization_structure<'a, T, H>(
210        &self,
211        add_args: &'a [T],
212        hedge_data: &H,
213        flow: Flow,
214    ) -> Option<Canonicalized<ShadowedStructure<Aind>>>
215    where
216        &'a T: Into<AtomOrView<'a>>,
217        H: HedgePolarizationData,
218    {
219        self.pol_symbol(flow).map(|name| {
220            ShadowedStructure::from_iter(
221                hedge_data.edge_spin_slots(),
222                name,
223                Some(add_args.iter().map(|arg| arg.into().into_owned()).collect()),
224            )
225        })
226    }
227}
228
229impl ReversibleEdge for EdgeData<&ParseEdge> {
230    fn pdg(&self) -> isize {
231        if let Some(p) = self.data.particle.particle() {
232            match self.orientation {
233                Orientation::Default | Orientation::Undirected => p.pdg_code,
234                Orientation::Reversed => p.pdg_code,
235            }
236        } else {
237            0
238        }
239    }
240
241    fn pol_symbol(&self, flow: Flow) -> Option<Symbol> {
242        if let Some(p) = self.data.particle.particle() {
243            // println!(
244            //     "pdg:{},orientation_aware:{},spin:{},flow:{:?}",
245            //     p.pdg_code,
246            //     self.pdg(),
247            //     p.spin,
248            //     flow
249            // );
250            match (p.spin, flow) {
251                (2, Flow::Sink) => {
252                    if self.pdg() > 0 {
253                        // println!("u");
254                        Some(GS.u)
255                    } else {
256                        // println!("vbar");
257                        Some(GS.vbar)
258                    }
259                }
260                (2, Flow::Source) => {
261                    if self.pdg() > 0 {
262                        // println!("ubar");
263                        Some(GS.ubar)
264                    } else {
265                        // println!("v");
266                        Some(GS.v)
267                    }
268                }
269                (3, Flow::Sink) => Some(GS.epsilon),
270                (3, Flow::Source) => Some(GS.epsilonbar),
271                _ => None,
272            }
273        } else {
274            None
275        }
276    }
277
278    fn polarization<'a, T, H>(&self, add_args: &'a [T], hedge_data: &H, flow: Flow) -> Atom
279    where
280        &'a T: Into<AtomOrView<'a>>,
281        H: HedgePolarizationData,
282    {
283        if let Some(name) = self.pol_symbol(flow) {
284            hedge_data.polarization(FunctionBuilder::new(name).add_args(add_args))
285        } else {
286            Atom::num(1)
287        }
288    }
289
290    fn polarization_structure<'a, T, H>(
291        &self,
292        add_args: &'a [T],
293        hedge_data: &H,
294        flow: Flow,
295    ) -> Option<Canonicalized<ShadowedStructure<Aind>>>
296    where
297        &'a T: Into<AtomOrView<'a>>,
298        H: HedgePolarizationData,
299    {
300        self.pol_symbol(flow).map(|name| {
301            ShadowedStructure::from_iter(
302                hedge_data.edge_spin_slots(),
303                name,
304                Some(add_args.iter().map(|arg| arg.into().into_owned()).collect()),
305            )
306        })
307    }
308}
309
310impl ReversibleEdge for EdgeData<(usize, isize)> {
311    fn pdg(&self) -> isize {
312        match self.orientation {
313            Orientation::Default | Orientation::Undirected => self.data.1,
314            Orientation::Reversed => -self.data.1,
315        }
316    }
317
318    fn pol_symbol(&self, flow: Flow) -> Option<Symbol> {
319        match (self.data.0, flow) {
320            (2, Flow::Sink) => {
321                if self.pdg() > 0 {
322                    Some(GS.u)
323                } else {
324                    Some(GS.vbar)
325                }
326            }
327            (2, Flow::Source) => {
328                if self.pdg() > 0 {
329                    Some(GS.ubar)
330                } else {
331                    Some(GS.v)
332                }
333            }
334            (3, Flow::Sink) => Some(GS.epsilon),
335            (3, Flow::Source) => Some(GS.epsilonbar),
336            _ => None,
337        }
338    }
339
340    fn polarization<'a, T, H>(&self, add_args: &'a [T], hedge_data: &H, flow: Flow) -> Atom
341    where
342        &'a T: Into<AtomOrView<'a>>,
343        H: HedgePolarizationData,
344    {
345        if let Some(name) = self.pol_symbol(flow) {
346            hedge_data.polarization(FunctionBuilder::new(name).add_args(add_args))
347        } else {
348            Atom::num(1)
349        }
350    }
351
352    fn polarization_structure<'a, T, H>(
353        &self,
354        add_args: &'a [T],
355        hedge_data: &H,
356        flow: Flow,
357    ) -> Option<Canonicalized<ShadowedStructure<Aind>>>
358    where
359        &'a T: Into<AtomOrView<'a>>,
360        H: HedgePolarizationData,
361    {
362        self.pol_symbol(flow).map(|name| {
363            ShadowedStructure::from_iter(
364                hedge_data.edge_spin_slots(),
365                name,
366                Some(add_args.iter().map(|arg| arg.into().into_owned()).collect()),
367            )
368        })
369    }
370}
371
372impl ReversibleEdge for EdgeData<ArcParticle> {
373    fn pdg(&self) -> isize {
374        match self.orientation {
375            Orientation::Default | Orientation::Undirected => self.data.pdg_code,
376            Orientation::Reversed => -self.data.pdg_code,
377        }
378    }
379
380    fn pol_symbol(&self, flow: Flow) -> Option<Symbol> {
381        match (self.data.spin, flow) {
382            (2, Flow::Sink) => {
383                if self.pdg() > 0 {
384                    Some(GS.u)
385                } else {
386                    Some(GS.vbar)
387                }
388            }
389            (2, Flow::Source) => {
390                if self.pdg() > 0 {
391                    Some(GS.ubar)
392                } else {
393                    Some(GS.v)
394                }
395            }
396            (3, Flow::Sink) => Some(GS.epsilon),
397            (3, Flow::Source) => Some(GS.epsilonbar),
398            _ => None,
399        }
400    }
401
402    fn polarization<'a, T, H>(&self, add_args: &'a [T], hedge_data: &H, flow: Flow) -> Atom
403    where
404        &'a T: Into<AtomOrView<'a>>,
405        H: HedgePolarizationData,
406    {
407        if let Some(name) = self.pol_symbol(flow) {
408            hedge_data.polarization(FunctionBuilder::new(name).add_args(add_args))
409        } else {
410            Atom::num(1)
411        }
412    }
413
414    fn polarization_structure<'a, T, H>(
415        &self,
416        add_args: &'a [T],
417        hedge_data: &H,
418        flow: Flow,
419    ) -> Option<Canonicalized<ShadowedStructure<Aind>>>
420    where
421        &'a T: Into<AtomOrView<'a>>,
422        H: HedgePolarizationData,
423    {
424        self.pol_symbol(flow).map(|name| {
425            ShadowedStructure::from_iter(
426                hedge_data.edge_spin_slots(),
427                name,
428                Some(add_args.iter().map(|arg| arg.into().into_owned()).collect()),
429            )
430        })
431    }
432}
433
434#[cfg(test)]
435mod test {
436
437    // use env_logger::WriteStyle;
438    use super::ReversibleEdge;
439    use idenso::{color::ColorSimplifier, dirac::GammaSimplifier, tensor::SymbolicNetParse};
440    use linnet::half_edge::{
441        involution::{EdgeData, Hedge},
442        subgraph::{SuBitGraph, SubSetLike},
443    };
444
445    use spenso::network::parsing::ParseSettings;
446    use symbolica::{
447        atom::{Atom, AtomCore},
448        parse_lit,
449    };
450
451    use crate::{
452        dot,
453        graph::{Graph, edge::ParseEdge, parse::IntoGraph},
454        initialisation::{initialise, test_initialise},
455        model::InputParamCard,
456        numerator::aind::Aind,
457        processes::{Amplitude, AmplitudeGraph, DotExportSettings},
458        settings::{
459            GlobalSettings, RuntimeSettings, global::GenerationSettings,
460            runtime::LockedRuntimeSettings,
461        },
462        utils::GS,
463        uv::UltravioletGraph,
464    };
465
466    #[test]
467    fn qqx_aaa_tree() {
468        test_initialise().unwrap();
469
470        let mut graph:AmplitudeGraph = dot!(digraph qqx_aaa_tree_1 {
471                    num="spenso::g(spenso::dind(spenso::cof(3, hedge(1))), spenso::cof(3, hedge(2)))/3"
472                    ext    [style=invis]
473                    ext -> v1:1 [particle="d" id=1];
474                    ext -> v3:2 [particle="d~" id=2];
475                    v1:3 -> ext [particle="a" id=3];
476                    v2:4 -> ext [particle="a" id=4];
477                    v3:0 -> ext [particle="a" id=0];
478                    v1 -> v2 [particle="d" id=5];
479                    v2 -> v3 [particle="d" id=6];
480        }).unwrap();
481
482        let vk = crate::utils::vakint().unwrap();
483
484        // let model = crate::utils::load_generic_model("sm");
485
486        graph.generate_cff(&GenerationSettings::default()).unwrap();
487        graph
488            .build_integrands(&GenerationSettings::default(), vk)
489            .unwrap();
490
491        println!("{}", graph.derived_data.all_mighty_integrand);
492    }
493
494    #[test]
495    fn simplify() {
496        test_initialise().unwrap();
497        let expr = parse_lit!(
498            -GC_1 * GC_11
499                ^ 2 * Q(4, mink(dim, gammalooprs::edge(4, 1)))
500                    * Q(5, mink(dim, gammalooprs::edge(5, 1)))
501                    * g(
502                        mink(dim, gammalooprs::hedge(16)),
503                        mink(dim, gammalooprs::hedge(17))
504                    )
505                    * g(
506                        coad(8, gammalooprs::hedge(16)),
507                        coad(8, gammalooprs::hedge(17))
508                    )
509                    * g(
510                        dind(cof(3, gammalooprs::hedge(4))),
511                        cof(3, gammalooprs::hedge(6))
512                    )
513                    * g(
514                        dind(cof(3, gammalooprs::hedge(5))),
515                        cof(3, gammalooprs::hedge(4))
516                    )
517                    * g(
518                        dind(cof(3, gammalooprs::hedge(6))),
519                        cof(3, gammalooprs::hedge(7))
520                    )
521                    * gamma(
522                        bis(4, gammalooprs::hedge(4)),
523                        bis(4, gammalooprs::hedge(6)),
524                        mink(dim, gammalooprs::hedge(8))
525                    )
526                    * gamma(
527                        bis(4, gammalooprs::hedge(7)),
528                        bis(4, gammalooprs::hedge(13)),
529                        mink(dim, gammalooprs::hedge(17))
530                    )
531                    * gamma(
532                        bis(4, gammalooprs::hedge(11)),
533                        bis(4, gammalooprs::hedge(5)),
534                        mink(dim, gammalooprs::hedge(16))
535                    )
536                    * gamma(
537                        bis(4, gammalooprs::hedge(5)),
538                        bis(4, gammalooprs::hedge(4)),
539                        mink(dim, gammalooprs::edge(4, 1))
540                    )
541                    * gamma(
542                        bis(4, gammalooprs::hedge(6)),
543                        bis(4, gammalooprs::hedge(7)),
544                        mink(dim, gammalooprs::edge(5, 1))
545                    )
546                    * t(
547                        coad(8, gammalooprs::hedge(16)),
548                        cof(3, gammalooprs::hedge(5)),
549                        dind(cof(3, gammalooprs::hedge(11)))
550                    )
551                    * t(
552                        coad(8, gammalooprs::hedge(17)),
553                        cof(3, gammalooprs::hedge(13)),
554                        dind(cof(3, gammalooprs::hedge(7)))
555                    ),
556            default_namespace = "spenso"
557        );
558
559        let net = expr
560            .parse_to_symbolic_net::<Aind>(&ParseSettings::default())
561            .unwrap();
562        println!("{}", net.dot_pretty());
563
564        let _ = expr.simplify_gamma();
565    }
566
567    #[test]
568    fn pols() {
569        initialise().unwrap();
570        let mut graphs: Vec<Graph> = dot!(
571            digraph bxatobx{
572                graph [
573                    num="v(0,spenso::bis(4,hedge(0)))*vbar(2,spenso::bis(4,hedge(2)))*ϵ(1,spenso::mink(4,hedge(1)))"
574                ]
575                ext    [style=invis]
576                ext -> A:1   [particle=a id=1]
577                ext -> A:2    [dir=back particle="b~" id=2]
578                A:0  -> ext  [dir=back particle="b~" id=0]
579            }
580
581            digraph batob{
582                graph [
583                    num="ubar(0,spenso::bis(4,hedge(0)))*u(2,spenso::bis(4,hedge(2)))*ϵ(1,spenso::mink(4,hedge(1)))"
584                ]
585                ext    [style=invis]
586                ext -> A:1   [particle=a id=1]
587                ext -> A:2    [particle="b" id=2]
588                A:0  -> ext  [particle="b" id=0]
589            }
590
591            digraph bbato{
592                graph [
593                    num="vbar(0,spenso::bis(4,hedge(0)))*u(2,spenso::bis(4,hedge(2)))*ϵ(1,spenso::mink(4,hedge(1)))"
594                ]
595                ext    [style=invis]
596                ext -> A:1   [particle=a id=1]
597                ext -> A:2    [particle="b" id=2]
598                ext -> A:0  [dir=back pdg=-5 id=0]
599            }
600
601            digraph bbato{
602                graph [
603                    num="vbar(0,spenso::bis(4,hedge(0)))*u(2,spenso::bis(4,hedge(2)))*ϵbar(1,spenso::mink(4,hedge(1)))"
604                ]
605                ext    [style=invis]
606                A:1 -> ext  [particle=a id=1]
607                ext -> A:2    [particle="b" id=2]
608                ext -> A:0  [dir=back pdg=-5 id=0]
609            }
610
611        )
612        .unwrap();
613
614        let model = crate::utils::load_generic_model("sm");
615        // Incoming (u, vbar) and outgoing (ubar, v) fermion states are fixed
616        // by the physical species and incidence, independently of arrow style.
617        let reference_symbols = [
618            [GS.v, GS.epsilon, GS.vbar],
619            [GS.ubar, GS.epsilon, GS.u],
620            [GS.vbar, GS.epsilon, GS.u],
621            [GS.vbar, GS.epsilonbar, GS.u],
622        ];
623        for (g, symbols) in graphs.iter_mut().zip(reference_symbols) {
624            let external: SuBitGraph = g.underlying.external_filter();
625            for hedge in external.included_iter() {
626                let eid = g.underlying[&hedge];
627                let edge = g.underlying.get_edge_data_full(hedge);
628                let parsed = ParseEdge::from(edge.data);
629                let parsed = EdgeData::new(&parsed, edge.orientation);
630                let flow = g.underlying.flow(hedge);
631                assert_eq!(parsed.pol_symbol(flow), Some(symbols[eid.0]));
632                assert_eq!(edge.pol_symbol(flow), parsed.pol_symbol(flow));
633            }
634            let pols = g.generate_polarizations();
635            assert_eq!(pols, g.global_prefactor.projector);
636            let serialized = g.dot_serialize(&DotExportSettings::default());
637            // This structural fixture stores its literal polarization reference
638            // as `num`; it does not define a complete scalar tensor network.
639            let round_trip = Graph::from_dot(
640                linnet::parser::DotGraph::from_string(serialized).unwrap(),
641                &model,
642            )
643            .unwrap();
644            assert_eq!(round_trip.generate_polarizations(), pols);
645            assert_eq!(round_trip.global_prefactor.projector, pols);
646            for hedge in external.included_iter() {
647                let original = g.underlying.get_edge_data_full(hedge);
648                let reparsed = round_trip.underlying.get_edge_data_full(hedge);
649                assert_eq!(
650                    reparsed.data.particle().map(|particle| particle.pdg_code),
651                    original.data.particle().map(|particle| particle.pdg_code),
652                );
653                assert_eq!(reparsed.orientation, original.orientation);
654            }
655            println!("{pols}");
656            let expected = &g.global_prefactor.num;
657            println!("{expected}");
658            assert_eq!(
659                expected, &pols,
660                "\nExpected: {:}\nActual: {:} for graph {:}",
661                expected, pols, g.name,
662            );
663        }
664    }
665
666    #[test]
667    fn user_numerator_tree_fixtures_match_ufo_tree_algebra() {
668        test_initialise().unwrap();
669        let mut model = crate::utils::load_generic_model("sm");
670        let mut parameters = InputParamCard::from_str(
671            include_str!("../../../../assets/models/json/sm/restrict_light_quarks_massless.json")
672                .to_owned(),
673            "json",
674        )
675        .unwrap();
676        model.simplify(&mut parameters).unwrap();
677        let mut expected = None;
678        for source in [
679            include_str!("../../../../tests/resources/graphs/qqx_aaa_tree_unprocessed.dot"),
680            include_str!(
681                "../../../../tests/resources/graphs/qqx_aaa_tree_user_numerator_unprocessed.dot"
682            ),
683            include_str!("../../../../tests/resources/graphs/qqx_aaa_tree_user_numerator.dot"),
684            include_str!("../../../../tests/resources/graphs/trees/qqx_aaa.dot"),
685        ] {
686            let graph: Graph = source.into_graph(&model).unwrap();
687            let numerator = (graph
688                .numerator(&graph.full_filter(), &graph.empty_subgraph())
689                .state
690                .expr
691                * &graph.global_prefactor.projector
692                * &graph.global_prefactor.num
693                * &graph.overall_factor)
694                .simplify_color();
695            if let Some(expected) = &expected {
696                assert_eq!(
697                    &numerator, expected,
698                    "manual tree must include both UFO propagator phases"
699                );
700            } else {
701                expected = Some(numerator);
702            }
703        }
704    }
705
706    #[test]
707    fn user_pentabox_numerator_closes_all_external_indices() {
708        test_initialise().unwrap();
709        let model = crate::utils::load_generic_model("sm");
710        let _: Graph =
711            include_str!("../../../../tests/resources/graphs/qqx_aaa_pentabox_user_numerator.dot")
712                .into_graph(&model)
713                .unwrap();
714    }
715
716    #[test]
717    fn tree() {
718        initialise().unwrap();
719        let mut model = crate::utils::load_generic_model("sm");
720        let mut parameters = InputParamCard::from_str(
721            include_str!("../../../../assets/models/json/sm/restrict_light_quarks_massless.json")
722                .to_owned(),
723            "json",
724        )
725        .unwrap();
726        model.simplify(&mut parameters).unwrap();
727        // The massless UFO reference includes i*pslash for each internal
728        // fermion propagator, matching the restricted SM Feynman rules.
729        let mut graphs: Vec<Graph> = dot!(
730            digraph qqx_aaa_tree_1 {
731                        num="spenso::g(spenso::dind(spenso::cof(3, hedge(1))), spenso::cof(3, hedge(2)))/3"
732                        ext    [style=invis]
733                        ext -> v1:1 [particle="d" id=1];
734                        ext -> v3:2 [dir=back particle="d~" id=2];
735                        v1:3 -> ext [particle="a" id=3];
736                        v2:4 -> ext [particle="a" id=4];
737                        v3:0 -> ext [particle="a" id=0];
738                        v1 -> v2 [particle="d" id=5];
739                        v2 -> v3 [particle="d" id=6];
740            }
741
742            digraph qqx_aaa_tree_1 {
743                        ext    [style=invis]
744                        ext -> v1:1 [particle="d" id=1];
745                        ext -> v3:2 [dir=back particle="d~" id=2];
746                        v1:3 -> ext [particle="a" id=3];
747                        v2:4 -> ext [particle="a" id=4];
748                        v3:0 -> ext [particle="a" id=0];
749                        v1 -> v2 [particle="d" id=5];
750                        v2 -> v3 [particle="d" id=6];
751                        num="spenso::g(spenso::dind(spenso::cof(3, hedge(1))), spenso::cof(3, hedge(2)))/3"
752            }
753
754
755            digraph qqx_aaa_tree_1_glob {
756            ext [style=invis];
757            ext -> v1:1 [particle="d", id=1];
758            ext -> v3:2 [dir=back particle="d~", id=2];
759            v1:3 -> ext [particle="a", id=3];
760            v2:4 -> ext [particle="a", id=4];
761            v3:0 -> ext [particle="a", id=0];
762            v1 -> v2 [particle="d", id=5];
763            v2 -> v3 [particle="d", id=6];
764            num=" UFO::GC_1^3
765                *spenso::g(spenso::cof(3,hedge(1)),spenso::dind(spenso::cof(3,hedge(5))))
766                *spenso::gamma(spenso::bis(4,hedge(5)),spenso::bis(4,hedge(1)),spenso::mink(4,hedge(3)))
767
768                *spenso::g(spenso::cof(3,hedge(5)),spenso::dind(spenso::cof(3,hedge(6))))
769                *1i*Q(5,spenso::mink(4,edge(5,1)))
770                *spenso::gamma(spenso::bis(4,hedge(6)),spenso::bis(4,hedge(5)),spenso::mink(4,edge(5,1)))
771
772
773                *spenso::g(spenso::cof(3,hedge(6)),spenso::dind(spenso::cof(3,hedge(7))))
774                *spenso::gamma(spenso::bis(4,hedge(7)),spenso::bis(4,hedge(6)),spenso::mink(4,hedge(4)))
775
776                *spenso::g(spenso::cof(3,hedge(7)),spenso::dind(spenso::cof(3,hedge(8))))
777                *spenso::gamma(spenso::bis(4,hedge(8)),spenso::bis(4,hedge(7)),spenso::mink(4,edge(6,1)))
778                *1i*Q(6,spenso::mink(4,edge(6,1)))
779
780                *spenso::g(spenso::cof(3,hedge(8)),spenso::dind(spenso::cof(3,hedge(2))))
781                *spenso::gamma(spenso::bis(4,hedge(2)),spenso::bis(4,hedge(8)),spenso::mink(4,hedge(0)))
782
783                   "
784             overall_factor="1"
785                projector="u(1,spenso::bis(4,hedge(1)))
786            *vbar(2,spenso::bis(4,hedge(2)))
787            *ϵbar(0,spenso::mink(4,hedge(0)))
788            *ϵbar(3,spenso::mink(4,hedge(3)))
789            *ϵbar(4,spenso::mink(4,hedge(4)))
790            *spenso::g(spenso::cof(3,hedge(2)),spenso::dind(spenso::cof(3,hedge(1))))/3"
791
792            edge [num="1"];
793            node [num="1"];
794            },
795            &model
796        )
797        .unwrap();
798
799        let mut a: Option<Atom> = None;
800        for g in &mut graphs {
801            let mut out = String::new();
802
803            let new_a = (g
804                .numerator(&g.full_filter(), &g.empty_subgraph())
805                .state
806                .expr
807                * &g.global_prefactor.projector
808                * &g.global_prefactor.num
809                * &g.overall_factor)
810                .simplify_color();
811            println!("New:{new_a}");
812            if let Some(a) = &a {
813                println!("Old:{a}");
814                assert_eq!(&new_a, a, "{}", (&new_a / a).expand().to_canonical_string());
815            } else {
816                a = Some(new_a);
817            }
818            g.dot_serialize_fmt(&mut out, &DotExportSettings::default())
819                .unwrap();
820            println!("{}", out);
821        }
822
823        let mut a = Amplitude::from_graph_list("test", graphs.clone()).unwrap();
824
825        let generation_pool = rayon::ThreadPoolBuilder::new()
826            .num_threads(1)
827            .build()
828            .unwrap();
829
830        let default_runtime_settings = RuntimeSettings::default();
831        let locked_runtime_settings = LockedRuntimeSettings::from(&default_runtime_settings);
832        a.preprocess(
833            &model,
834            &GenerationSettings::default(),
835            &locked_runtime_settings,
836            &generation_pool,
837        )
838        .unwrap();
839    }
840
841    #[test]
842    fn evaluate_pols() {
843        initialise().unwrap();
844        let model = crate::utils::load_generic_model("sm");
845
846        let mut graph = Graph::from_dot(
847            linnet::dot!(
848                digraph bxatobx {
849                    bla [style=invis]
850                    bla -> A:1 [particle=a id=0]
851                    bla -> A:2 [particle="b~" id=2]
852                    A:0 -> bla [particle="b~" id=1]
853                }
854            )
855            .unwrap(),
856            &model,
857        )
858        .unwrap();
859        graph.global_prefactor.num *=
860            graph.underlying[Hedge(0)].color_kronekers(&graph.underlying[Hedge(2)]);
861        graph.validate_full_numerator_tensor_network().unwrap();
862        let mut amp = Amplitude::from_graph_list("name", vec![graph]).unwrap();
863
864        let proc_set = GenerationSettings::default();
865        let thread_pool = rayon::ThreadPoolBuilder::new()
866            .num_threads(1)
867            .build()
868            .unwrap();
869
870        let default_runtime_settings = RuntimeSettings::default();
871        let locked_runtime_settings = LockedRuntimeSettings::from(&default_runtime_settings);
872        amp.preprocess(&model, &proc_set, &locked_runtime_settings, &thread_pool)
873            .unwrap();
874
875        amp.build_integrand(
876            &model,
877            "test_process",
878            &GlobalSettings::default(),
879            (&RuntimeSettings::default()).into(),
880            &thread_pool,
881        )
882        .unwrap();
883
884        // integrand.evaluate_sample(sample, wgt, iter, use_f128, max_eval);
885    }
886
887    mod failing {
888        use super::*;
889
890        #[test]
891        fn vertex_rule() {
892            let mut graphs: Vec<Graph> = dot!(
893                digraph dxda{
894                    ext [style=invis]
895                    ext->v1:0[particle="d" id=0]
896                    ext->v1:1[particle="d~" id=1]
897                    v1:2->ext[particle="a" id=2]
898                }
899
900            )
901            .unwrap();
902
903            for g in &mut graphs {
904                let mut out = String::new();
905                g.dot_serialize_fmt(&mut out, &DotExportSettings::default())
906                    .unwrap();
907                println!("{}", out);
908
909                assert!(g.iter_nodes().all(|(_, _, v)| {
910                    let a = v.vertex_rule.as_ref().unwrap().name.as_str();
911                    a == "V_71"
912                }),);
913            }
914        }
915
916        #[test]
917        fn pslash() {
918            let mut graphs: Vec<Graph> = dot!(
919                digraph dxda{
920                    ext [style=invis]
921                    node[num=1]
922                    ext->v1:0[particle="d" id=0]
923                    v1:1->v2:2[particle="d" id=1]
924                    v2:3->ext[particle="d" id=2]
925                }
926
927            )
928            .unwrap();
929
930            for g in &mut graphs {
931                let mut out = String::new();
932                g.dot_serialize_fmt(&mut out, &DotExportSettings::default())
933                    .unwrap();
934                println!("{}", out);
935            }
936        }
937
938        #[test]
939        fn dod_override() {
940            let _gr:Vec<Graph> = dot!(digraph g{
941                node [num=1]
942                edge [num=1]
943                a->b[dod=-100]
944                b->c[dod="-100"]
945            }
946            digraph triangle_ct {
947            num="(-1*_gammaloop::P(1,spenso::cind(1))*_gammaloop::P(2,spenso::cind(1))+-1*_gammaloop::P(1,spenso::cind(2))*_gammaloop::P(2,spenso::cind(2))+-1*_gammaloop::P(1,spenso::cind(3))*_gammaloop::P(2,spenso::cind(3))+_gammaloop::P(1,spenso::cind(0))*_gammaloop::P(2,spenso::cind(0)))^-2*(-1*_gammaloop::P(2,spenso::mink(4,python::mu2))+-1*_gammaloop::Q(6,spenso::mink(4,python::mu2))+_gammaloop::P(1,spenso::mink(4,python::mu7)))*(-1*_gammaloop::Q(7,spenso::mink(4,python::mu2))+_gammaloop::P(2,spenso::mink(4,python::mu2)))*-1/4*_gammaloop::G^2*_gammaloop::P(1,spenso::mink(4,python::mu3))*_gammaloop::P(1,spenso::mink(4,python::mu5))*_gammaloop::P(2,spenso::mink(4,python::mu4))*_gammaloop::P(2,spenso::mink(4,python::mu6))*spenso::gamma(spenso::bis(4,_gammaloop::hedge(2)),spenso::bis(4,python::s1),spenso::mink(4,python::mu1))*spenso::gamma(spenso::bis(4,python::s1),spenso::bis(4,python::s2),spenso::mink(4,python::mu2))*spenso::gamma(spenso::bis(4,python::s2),spenso::bis(4,python::s3),spenso::mink(4,python::mu3))*spenso::gamma(spenso::bis(4,python::s3),spenso::bis(4,python::tree_form_factor_spinor_2),spenso::mink(4,python::mu4))*spenso::gamma(spenso::bis(4,python::s4),spenso::bis(4,python::s6),spenso::mink(4,python::mu7))*spenso::gamma(spenso::bis(4,python::s5),spenso::bis(4,python::s4),spenso::mink(4,python::mu6))*spenso::gamma(spenso::bis(4,python::s6),spenso::bis(4,_gammaloop::hedge(1)),spenso::mink(4,python::mu1))*spenso::gamma(spenso::bis(4,python::tree_form_factor_spinor_1),spenso::bis(4,python::s5),spenso::mink(4,python::mu5))";
948            overall_factor="1";
949            projector="((-1*_gammaloop::P(3,spenso::cind(0))+-1*_gammaloop::P(4,spenso::cind(0))+_gammaloop::P(1,spenso::cind(0)))^2+(-1*_gammaloop::P(3,spenso::cind(1))+-1*_gammaloop::P(4,spenso::cind(1))+_gammaloop::P(1,spenso::cind(1)))^2*-1+(-1*_gammaloop::P(3,spenso::cind(2))+-1*_gammaloop::P(4,spenso::cind(2))+_gammaloop::P(1,spenso::cind(2)))^2*-1+(-1*_gammaloop::P(3,spenso::cind(3))+-1*_gammaloop::P(4,spenso::cind(3))+_gammaloop::P(1,spenso::cind(3)))^2*-1)^-1*((-1*_gammaloop::P(3,spenso::cind(0))+_gammaloop::P(1,spenso::cind(0)))^2+(-1*_gammaloop::P(3,spenso::cind(1))+_gammaloop::P(1,spenso::cind(1)))^2*-1+(-1*_gammaloop::P(3,spenso::cind(2))+_gammaloop::P(1,spenso::cind(2)))^2*-1+(-1*_gammaloop::P(3,spenso::cind(3))+_gammaloop::P(1,spenso::cind(3)))^2*-1)^-1*(-1*_gammaloop::P(3,spenso::mink(4,_gammaloop::edge(5,1)))+_gammaloop::P(1,spenso::mink(4,_gammaloop::edge(5,1))))*(-1*_gammaloop::P(3,spenso::mink(4,_gammaloop::edge(6,1)))+-1*_gammaloop::P(4,spenso::mink(4,_gammaloop::edge(6,1)))+_gammaloop::P(1,spenso::mink(4,_gammaloop::edge(6,1))))*-1/27*_gammaloop::ee^3*_gammaloop::u(1,spenso::bis(4,_gammaloop::hedge(1)))*_gammaloop::vbar(2,spenso::bis(4,_gammaloop::hedge(2)))*_gammaloop::ϵbar(0,spenso::mink(4,_gammaloop::hedge(0)))*_gammaloop::ϵbar(3,spenso::mink(4,_gammaloop::hedge(3)))*_gammaloop::ϵbar(4,spenso::mink(4,_gammaloop::hedge(4)))*spenso::gamma(spenso::bis(4,_gammaloop::hedge(5)),spenso::bis(4,python::tree_form_factor_spinor_1),spenso::mink(4,_gammaloop::hedge(3)))*spenso::gamma(spenso::bis(4,_gammaloop::hedge(6)),spenso::bis(4,_gammaloop::hedge(5)),spenso::mink(4,_gammaloop::edge(5,1)))*spenso::gamma(spenso::bis(4,_gammaloop::hedge(7)),spenso::bis(4,_gammaloop::hedge(6)),spenso::mink(4,_gammaloop::hedge(4)))*spenso::gamma(spenso::bis(4,_gammaloop::hedge(8)),spenso::bis(4,_gammaloop::hedge(7)),spenso::mink(4,_gammaloop::edge(6,1)))*spenso::gamma(spenso::bis(4,python::tree_form_factor_spinor_2),spenso::bis(4,_gammaloop::hedge(8)),spenso::mink(4,_gammaloop::hedge(0)))";
950            edge [num="1", dod="-100"];
951            node [num="1", dod="-100"];
952            ext [style=invis];
953            ext -> vl1:1 [particle="d", id=1];
954            ext -> vl2:2 [dir=back particle="d~", id=2];
955            v1:3 -> ext [id=3];
956            v1:4 -> ext [id=4];
957            v1:0 -> ext [id=0];
958            vl1 -> v1 [particle="d", id=5];
959            v1 -> vl2 [particle="d", id=6];
960            vl1 -> vl2 [particle="g", id=7, lmb_id=0];
961            }
962    )
963            .unwrap();
964        }
965    }
966}
gammalooprs/numerator/graph.rs

gammalooprs/numerator/
graph.rs
