spenso_hep_lib/
lib.rs

1use std::{ops::Neg, sync::LazyLock};
2
3use idenso::{color::CS, dirac::AGS, representations::initialize};
4
5use spenso::{
6    algebra::complex::Complex,
7    network::{
8        Network,
9        library::{
10            TensorLibraryData,
11            function_lib::{INBUILTS, PanicMissingConcrete, SymbolLib},
12            symbolic::{ExplicitKey, TensorLibrary},
13        },
14        parsing::ShadowedStructure,
15        store::NetworkStore,
16    },
17    structure::{PermutedStructure, TensorStructure, abstract_index::AbstractIndex, slot::AbsInd},
18    tensors::{
19        complex::RealOrComplexTensor,
20        data::{SetTensorData, SparseTensor, StorageTensor},
21        parametric::{MixedTensor, ParamOrConcrete},
22    },
23};
24use symbolica::{
25    atom::{Atom, Symbol},
26    parse_lit,
27};
28
29#[allow(clippy::similar_names)]
30pub fn gamma_data_dirac<T, N>(structure: N, one: T, zero: T) -> SparseTensor<Complex<T>, N>
31where
32    T: Clone + Neg<Output = T>,
33    N: TensorStructure,
34{
35    let c1 = Complex::<T>::new(one.clone(), zero.clone());
36    let z = Complex::<T>::new(zero.clone(), zero.clone());
37    let cn1 = Complex::<T>::new(-one.clone(), zero.clone());
38    let ci = Complex::<T>::new(zero.clone(), one.clone());
39    let cni = Complex::<T>::new(zero.clone(), -one.clone());
40    let mut gamma = SparseTensor::empty(structure, z);
41    // ! No check on actual structure, should expext mink,bis,bis
42
43    // dirac gamma matrices
44
45    gamma.set(&[0, 0, 0], c1.clone()).unwrap();
46    gamma.set(&[0, 1, 1], c1.clone()).unwrap();
47    gamma.set(&[0, 2, 2], cn1.clone()).unwrap();
48    gamma.set(&[0, 3, 3], cn1.clone()).unwrap();
49
50    gamma.set(&[1, 0, 3], c1.clone()).unwrap();
51    gamma.set(&[1, 1, 2], c1.clone()).unwrap();
52    gamma.set(&[1, 2, 1], cn1.clone()).unwrap();
53    gamma.set(&[1, 3, 0], cn1.clone()).unwrap();
54
55    gamma.set(&[2, 0, 3], cni.clone()).unwrap();
56    gamma.set(&[2, 1, 2], ci.clone()).unwrap();
57    gamma.set(&[2, 2, 1], ci.clone()).unwrap();
58    gamma.set(&[2, 3, 0], cni.clone()).unwrap();
59
60    gamma.set(&[3, 0, 2], c1.clone()).unwrap();
61    gamma.set(&[3, 1, 3], cn1.clone()).unwrap();
62    gamma.set(&[3, 2, 0], cn1.clone()).unwrap();
63    gamma.set(&[3, 3, 1], c1.clone()).unwrap();
64
65    gamma //.to_dense()
66}
67
68#[allow(clippy::similar_names)]
69pub fn gamma_data_weyl<T, N>(structure: N, one: T, zero: T) -> SparseTensor<Complex<T>, N>
70where
71    T: Neg<Output = T> + Clone,
72    N: TensorStructure,
73{
74    let z = Complex::<T>::new(zero.clone(), zero.clone());
75    let c1 = Complex::<T>::new(one.clone(), zero.clone());
76    let cn1 = Complex::<T>::new(-one.clone(), zero.clone());
77    let ci = Complex::<T>::new(zero.clone(), one.clone());
78    let cni = Complex::<T>::new(zero.clone(), -one.clone());
79    let mut gamma = SparseTensor::empty(structure, z);
80    // ! No check on actual structure, should expext mink,bis,bis
81
82    // dirac gamma matrices
83
84    gamma.set(&[0, 2, 0], c1.clone()).unwrap();
85    gamma.set(&[1, 3, 0], c1.clone()).unwrap();
86    gamma.set(&[2, 0, 0], c1.clone()).unwrap();
87    gamma.set(&[3, 1, 0], c1.clone()).unwrap();
88
89    gamma.set(&[0, 3, 1], c1.clone()).unwrap();
90    gamma.set(&[1, 2, 1], c1.clone()).unwrap();
91    gamma.set(&[2, 1, 1], cn1.clone()).unwrap();
92    gamma.set(&[3, 0, 1], cn1.clone()).unwrap();
93
94    gamma.set(&[0, 3, 2], cni.clone()).unwrap();
95    gamma.set(&[1, 2, 2], ci.clone()).unwrap();
96    gamma.set(&[2, 1, 2], ci.clone()).unwrap();
97    gamma.set(&[3, 0, 2], cni.clone()).unwrap();
98
99    gamma.set(&[0, 2, 3], c1.clone()).unwrap();
100    gamma.set(&[1, 3, 3], cn1.clone()).unwrap();
101    gamma.set(&[2, 0, 3], cn1.clone()).unwrap();
102    gamma.set(&[3, 1, 3], c1.clone()).unwrap();
103
104    gamma //.to_dense()
105}
106
107#[allow(clippy::similar_names)]
108pub fn gamma_transpose_weyl<T, N>(structure: N, one: T, zero: T) -> SparseTensor<Complex<T>, N>
109where
110    T: Neg<Output = T> + Clone,
111    N: TensorStructure,
112{
113    let z = Complex::<T>::new(zero.clone(), zero.clone());
114    let c1 = Complex::<T>::new(one.clone(), zero.clone());
115    let cn1 = Complex::<T>::new(-one.clone(), zero.clone());
116    let ci = Complex::<T>::new(zero.clone(), one.clone());
117    let cni = Complex::<T>::new(zero.clone(), -one.clone());
118    let mut gamma = SparseTensor::empty(structure, z);
119    // ! No check on actual structure, should expext mink,bis,bis
120
121    // dirac gamma matrices
122
123    gamma.set(&[2, 0, 0], c1.clone()).unwrap();
124    gamma.set(&[3, 1, 0], c1.clone()).unwrap();
125    gamma.set(&[0, 2, 0], c1.clone()).unwrap();
126    gamma.set(&[1, 3, 0], c1.clone()).unwrap();
127
128    gamma.set(&[3, 0, 1], c1.clone()).unwrap();
129    gamma.set(&[2, 1, 1], c1.clone()).unwrap();
130    gamma.set(&[1, 2, 1], cn1.clone()).unwrap();
131    gamma.set(&[0, 3, 1], cn1.clone()).unwrap();
132
133    gamma.set(&[3, 0, 2], cni.clone()).unwrap();
134    gamma.set(&[2, 1, 2], ci.clone()).unwrap();
135    gamma.set(&[1, 2, 2], ci.clone()).unwrap();
136    gamma.set(&[0, 3, 2], cni.clone()).unwrap();
137
138    gamma.set(&[2, 0, 3], c1.clone()).unwrap();
139    gamma.set(&[3, 1, 3], cn1.clone()).unwrap();
140    gamma.set(&[0, 2, 3], cn1.clone()).unwrap();
141    gamma.set(&[1, 3, 3], c1.clone()).unwrap();
142
143    gamma //.to_dense()
144}
145
146#[allow(clippy::similar_names)]
147pub fn gamma_conj_data_weyl<T, N>(structure: N, one: T, zero: T) -> SparseTensor<Complex<T>, N>
148where
149    T: Neg<Output = T> + Clone,
150    N: TensorStructure,
151{
152    gamma_data_weyl(structure, one, zero).map_data(|a| {
153        let Complex { re, im } = a;
154        Complex { re, im: -im }
155    })
156}
157
158#[allow(clippy::similar_names)]
159pub fn gamma_adj_data_weyl<T, N>(structure: N, one: T, zero: T) -> SparseTensor<Complex<T>, N>
160where
161    T: Neg<Output = T> + Clone,
162    N: TensorStructure,
163{
164    gamma_transpose_weyl(structure, one, zero).map_data(|a| {
165        let Complex { re, im } = a;
166        Complex { re, im: -im }
167    })
168}
169
170pub fn gamma0_weyl<T, N>(structure: N, one: T, zero: T) -> SparseTensor<Complex<T>, N>
171where
172    T: Clone,
173    N: TensorStructure,
174{
175    let c1 = Complex::<T>::new(one, zero.clone());
176    let z = Complex::<T>::new(zero.clone(), zero.clone());
177    let mut gamma0 = SparseTensor::empty(structure, z);
178    // ! No check on actual structure, should expext bis,bis,lor
179
180    // dirac gamma0 matrices
181
182    gamma0.set(&[0, 2], c1.clone()).unwrap();
183    gamma0.set(&[1, 3], c1.clone()).unwrap();
184    gamma0.set(&[2, 0], c1.clone()).unwrap();
185    gamma0.set(&[3, 1], c1.clone()).unwrap();
186
187    gamma0
188}
189
190pub fn gamma5_dirac_data<T, N>(structure: N, one: T, zero: T) -> SparseTensor<Complex<T>, N>
191where
192    T: Clone,
193    N: TensorStructure,
194{
195    let c1 = Complex::<T>::new(one, zero.clone());
196
197    let z = Complex::<T>::new(zero.clone(), zero.clone());
198    let mut gamma5 = SparseTensor::empty(structure, z);
199
200    gamma5.set(&[0, 2], c1.clone()).unwrap();
201    gamma5.set(&[1, 3], c1.clone()).unwrap();
202    gamma5.set(&[2, 0], c1.clone()).unwrap();
203    gamma5.set(&[3, 1], c1.clone()).unwrap();
204
205    gamma5
206}
207
208pub fn gamma5_weyl_data<T, N>(structure: N, one: T, zero: T) -> SparseTensor<Complex<T>, N>
209where
210    T: Clone + Neg<Output = T>,
211    N: TensorStructure,
212{
213    let z = Complex::<T>::new(zero.clone(), zero.clone());
214    let c1 = Complex::<T>::new(one, zero);
215
216    let mut gamma5 = SparseTensor::empty(structure, z);
217
218    gamma5.set(&[0, 0], -c1.clone()).unwrap();
219    gamma5.set(&[1, 1], -c1.clone()).unwrap();
220    gamma5.set(&[2, 2], c1.clone()).unwrap();
221    gamma5.set(&[3, 3], c1.clone()).unwrap();
222
223    gamma5
224}
225
226#[allow(clippy::similar_names)]
227pub fn proj_m_data_dirac<T, N>(structure: N, half: T, zero: T) -> SparseTensor<Complex<T>, N>
228where
229    T: Clone + Neg<Output = T>,
230    N: TensorStructure,
231{
232    let z = Complex::<T>::new(zero.clone(), zero.clone());
233    // ProjM(1,2) Left chirality projector (( 1−γ5)/ 2 )_s1_s2
234
235    let chalf = Complex::<T>::new(half.clone(), zero.clone());
236    let cnhalf = Complex::<T>::new(-half, zero);
237
238    let mut proj_m = SparseTensor::empty(structure, z);
239
240    proj_m.set(&[0, 0], chalf.clone()).unwrap();
241    proj_m.set(&[1, 1], chalf.clone()).unwrap();
242    proj_m.set(&[2, 2], chalf.clone()).unwrap();
243    proj_m.set(&[3, 3], chalf.clone()).unwrap();
244
245    proj_m.set(&[0, 2], cnhalf.clone()).unwrap();
246    proj_m.set(&[1, 3], cnhalf.clone()).unwrap();
247    proj_m.set(&[2, 0], cnhalf.clone()).unwrap();
248    proj_m.set(&[3, 1], cnhalf.clone()).unwrap();
249
250    proj_m
251}
252
253#[allow(clippy::similar_names)]
254pub fn proj_m_data_weyl<T, N>(structure: N, one: T, zero: T) -> SparseTensor<Complex<T>, N>
255where
256    T: Clone,
257    N: TensorStructure,
258{
259    let z = Complex::<T>::new(zero.clone(), zero.clone());
260    // ProjM(1,2) Left chirality projector (( 1−γ5)/ 2 )_s1_s2
261    let c1 = Complex::<T>::new(one, zero);
262    let mut proj_m = SparseTensor::empty(structure, z);
263
264    proj_m.set(&[0, 0], c1.clone()).unwrap();
265    proj_m.set(&[1, 1], c1.clone()).unwrap();
266
267    proj_m
268}
269
270pub fn proj_p_data_dirac<T, N>(structure: N, half: T, zero: T) -> SparseTensor<Complex<T>, N>
271where
272    T: Clone,
273    N: TensorStructure,
274{
275    let z = Complex::<T>::new(zero.clone(), zero.clone());
276    // ProjP(1,2) Right chirality projector (( 1+γ5)/ 2 )_s1_s2
277    let chalf = Complex::<T>::new(half, zero);
278
279    let mut proj_p = SparseTensor::empty(structure, z);
280
281    proj_p
282        .set(&[0, 0], chalf.clone())
283        .unwrap_or_else(|_| unreachable!());
284    proj_p
285        .set(&[1, 1], chalf.clone())
286        .unwrap_or_else(|_| unreachable!());
287    proj_p
288        .set(&[2, 2], chalf.clone())
289        .unwrap_or_else(|_| unreachable!());
290    proj_p
291        .set(&[3, 3], chalf.clone())
292        .unwrap_or_else(|_| unreachable!());
293
294    proj_p
295        .set(&[0, 2], chalf.clone())
296        .unwrap_or_else(|_| unreachable!());
297    proj_p
298        .set(&[1, 3], chalf.clone())
299        .unwrap_or_else(|_| unreachable!());
300    proj_p
301        .set(&[2, 0], chalf.clone())
302        .unwrap_or_else(|_| unreachable!());
303    proj_p
304        .set(&[3, 1], chalf.clone())
305        .unwrap_or_else(|_| unreachable!());
306
307    proj_p
308}
309
310#[allow(clippy::similar_names)]
311pub fn proj_p_data_weyl<T, N>(structure: N, one: T, zero: T) -> SparseTensor<Complex<T>, N>
312where
313    T: Clone,
314    N: TensorStructure,
315{
316    let z = Complex::<T>::new(zero.clone(), zero.clone());
317    // ProjM(1,2) Left chirality projector (( 1−γ5)/ 2 )_s1_s2
318    let c1 = Complex::<T>::new(one, zero);
319    let mut proj_p = SparseTensor::empty(structure, z);
320
321    proj_p.set(&[2, 2], c1.clone()).unwrap();
322    proj_p.set(&[3, 3], c1.clone()).unwrap();
323
324    proj_p
325}
326
327/// Fundamental SU(3) generators in the normalization `Tr(T^a T^b)=1/2 delta^{ab}`.
328///
329/// The index order follows `CS.t_strct`: adjoint, fundamental, anti-fundamental.
330pub fn su3_generator_data<N>(structure: N) -> SparseTensor<Complex<f64>, N>
331where
332    N: TensorStructure,
333{
334    let z = Complex::new(0., 0.);
335    let mut t = SparseTensor::empty(structure, z);
336    let h = 0.5;
337    let s = 1. / (2. * 3_f64.sqrt());
338
339    t.set(&[0, 0, 1], Complex::new(h, 0.)).unwrap();
340    t.set(&[0, 1, 0], Complex::new(h, 0.)).unwrap();
341
342    t.set(&[1, 0, 1], Complex::new(0., -h)).unwrap();
343    t.set(&[1, 1, 0], Complex::new(0., h)).unwrap();
344
345    t.set(&[2, 0, 0], Complex::new(h, 0.)).unwrap();
346    t.set(&[2, 1, 1], Complex::new(-h, 0.)).unwrap();
347
348    t.set(&[3, 0, 2], Complex::new(h, 0.)).unwrap();
349    t.set(&[3, 2, 0], Complex::new(h, 0.)).unwrap();
350
351    t.set(&[4, 0, 2], Complex::new(0., -h)).unwrap();
352    t.set(&[4, 2, 0], Complex::new(0., h)).unwrap();
353
354    t.set(&[5, 1, 2], Complex::new(h, 0.)).unwrap();
355    t.set(&[5, 2, 1], Complex::new(h, 0.)).unwrap();
356
357    t.set(&[6, 1, 2], Complex::new(0., -h)).unwrap();
358    t.set(&[6, 2, 1], Complex::new(0., h)).unwrap();
359
360    t.set(&[7, 0, 0], Complex::new(s, 0.)).unwrap();
361    t.set(&[7, 1, 1], Complex::new(s, 0.)).unwrap();
362    t.set(&[7, 2, 2], Complex::new(-2. * s, 0.)).unwrap();
363
364    t
365}
366
367/// SU(3) structure constants for `[T^a,T^b]=i f^{abc} T^c`.
368pub fn su3_structure_f_data<N>(structure: N) -> SparseTensor<f64, N>
369where
370    N: TensorStructure,
371{
372    let mut f = SparseTensor::empty(structure, 0.);
373
374    fn set_antisymmetric<N>(f: &mut SparseTensor<f64, N>, a: usize, b: usize, c: usize, value: f64)
375    where
376        N: TensorStructure,
377    {
378        f.set(&[a, b, c], value).unwrap();
379        f.set(&[b, c, a], value).unwrap();
380        f.set(&[c, a, b], value).unwrap();
381        f.set(&[b, a, c], -value).unwrap();
382        f.set(&[a, c, b], -value).unwrap();
383        f.set(&[c, b, a], -value).unwrap();
384    }
385
386    set_antisymmetric(&mut f, 0, 1, 2, 1.);
387    set_antisymmetric(&mut f, 0, 3, 6, 0.5);
388    set_antisymmetric(&mut f, 0, 4, 5, -0.5);
389    set_antisymmetric(&mut f, 1, 3, 5, 0.5);
390    set_antisymmetric(&mut f, 1, 4, 6, 0.5);
391    set_antisymmetric(&mut f, 2, 3, 4, 0.5);
392    set_antisymmetric(&mut f, 2, 5, 6, -0.5);
393    set_antisymmetric(&mut f, 3, 4, 7, 3_f64.sqrt() / 2.);
394    set_antisymmetric(&mut f, 5, 6, 7, 3_f64.sqrt() / 2.);
395
396    f
397}
398
399/// Exact Atom-backed SU(3) generators, useful for symbolic tensor libraries.
400pub fn su3_generator_data_atom<N>(structure: N) -> SparseTensor<Atom, N>
401where
402    N: TensorStructure,
403{
404    let mut t = SparseTensor::empty(structure, Atom::Zero);
405    let half = Atom::num(1) / Atom::num(2);
406    let sqrt3 = parse_lit!(sqrt(3));
407    let t8 = sqrt3.clone() / Atom::num(6);
408
409    t.set(&[0, 0, 1], half.clone()).unwrap();
410    t.set(&[0, 1, 0], half.clone()).unwrap();
411
412    t.set(&[1, 0, 1], -Atom::i() * half.clone()).unwrap();
413    t.set(&[1, 1, 0], Atom::i() * half.clone()).unwrap();
414
415    t.set(&[2, 0, 0], half.clone()).unwrap();
416    t.set(&[2, 1, 1], -half.clone()).unwrap();
417
418    t.set(&[3, 0, 2], half.clone()).unwrap();
419    t.set(&[3, 2, 0], half.clone()).unwrap();
420
421    t.set(&[4, 0, 2], -Atom::i() * half.clone()).unwrap();
422    t.set(&[4, 2, 0], Atom::i() * half.clone()).unwrap();
423
424    t.set(&[5, 1, 2], half.clone()).unwrap();
425    t.set(&[5, 2, 1], half.clone()).unwrap();
426
427    t.set(&[6, 1, 2], -Atom::i() * half).unwrap();
428    t.set(&[6, 2, 1], Atom::i() / Atom::num(2)).unwrap();
429
430    t.set(&[7, 0, 0], t8.clone()).unwrap();
431    t.set(&[7, 1, 1], t8.clone()).unwrap();
432    t.set(&[7, 2, 2], -sqrt3 / Atom::num(3)).unwrap();
433
434    t
435}
436
437/// Exact Atom-backed SU(3) structure constants.
438pub fn su3_structure_f_data_atom<N>(structure: N) -> SparseTensor<Atom, N>
439where
440    N: TensorStructure,
441{
442    let mut f = SparseTensor::empty(structure, Atom::Zero);
443
444    fn set_antisymmetric<N>(
445        f: &mut SparseTensor<Atom, N>,
446        a: usize,
447        b: usize,
448        c: usize,
449        value: Atom,
450    ) where
451        N: TensorStructure,
452    {
453        f.set(&[a, b, c], value.clone()).unwrap();
454        f.set(&[b, c, a], value.clone()).unwrap();
455        f.set(&[c, a, b], value.clone()).unwrap();
456        f.set(&[b, a, c], -value.clone()).unwrap();
457        f.set(&[a, c, b], -value.clone()).unwrap();
458        f.set(&[c, b, a], -value).unwrap();
459    }
460
461    let half = Atom::num(1) / Atom::num(2);
462    let sqrt3_half = parse_lit!(sqrt(3)) / Atom::num(2);
463
464    set_antisymmetric(&mut f, 0, 1, 2, Atom::num(1));
465    set_antisymmetric(&mut f, 0, 3, 6, half.clone());
466    set_antisymmetric(&mut f, 0, 4, 5, -half.clone());
467    set_antisymmetric(&mut f, 1, 3, 5, half.clone());
468    set_antisymmetric(&mut f, 1, 4, 6, half.clone());
469    set_antisymmetric(&mut f, 2, 3, 4, half.clone());
470    set_antisymmetric(&mut f, 2, 5, 6, -half);
471    set_antisymmetric(&mut f, 3, 4, 7, sqrt3_half.clone());
472    set_antisymmetric(&mut f, 5, 6, 7, sqrt3_half);
473
474    f
475}
476
477#[allow(clippy::similar_names)]
478pub fn sigma_data<T, N>(structure: N, one: T, zero: T) -> SparseTensor<Complex<T>, N>
479where
480    T: Clone + Neg<Output = T>,
481    N: TensorStructure,
482{
483    let z = Complex::<T>::new(zero.clone(), zero.clone());
484    let c1 = Complex::<T>::new(one.clone(), zero.clone());
485    let cn1 = Complex::<T>::new(-one.clone(), zero.clone());
486    let ci = Complex::<T>::new(zero.clone(), one.clone());
487    let cni = Complex::<T>::new(zero.clone(), -one.clone());
488
489    let mut sigma = SparseTensor::empty(structure, z);
490    sigma.set(&[0, 2, 0, 1], c1.clone()).unwrap();
491    sigma.set(&[0, 2, 3, 0], c1.clone()).unwrap();
492    sigma.set(&[0, 3, 1, 2], c1.clone()).unwrap();
493    sigma.set(&[1, 0, 2, 2], c1.clone()).unwrap();
494    sigma.set(&[1, 1, 1, 2], c1.clone()).unwrap();
495    sigma.set(&[1, 3, 0, 2], c1.clone()).unwrap();
496    sigma.set(&[2, 2, 1, 0], c1.clone()).unwrap();
497    sigma.set(&[2, 2, 2, 1], c1.clone()).unwrap();
498    sigma.set(&[2, 3, 3, 2], c1.clone()).unwrap();
499    sigma.set(&[3, 0, 0, 2], c1.clone()).unwrap();
500    sigma.set(&[3, 3, 2, 2], c1.clone()).unwrap();
501    sigma.set(&[3, 1, 3, 2], c1.clone()).unwrap();
502    sigma.set(&[0, 1, 3, 0], ci.clone()).unwrap();
503    sigma.set(&[0, 3, 1, 1], ci.clone()).unwrap();
504    sigma.set(&[0, 3, 2, 0], ci.clone()).unwrap();
505    sigma.set(&[1, 0, 3, 3], ci.clone()).unwrap();
506    sigma.set(&[1, 1, 0, 3], ci.clone()).unwrap();
507    sigma.set(&[1, 1, 2, 0], ci.clone()).unwrap();
508    sigma.set(&[2, 1, 1, 0], ci.clone()).unwrap();
509    sigma.set(&[2, 3, 0, 0], ci.clone()).unwrap();
510    sigma.set(&[2, 3, 3, 1], ci.clone()).unwrap();
511    sigma.set(&[3, 0, 1, 3], ci.clone()).unwrap();
512    sigma.set(&[3, 1, 0, 0], ci.clone()).unwrap();
513    sigma.set(&[3, 1, 2, 3], ci.clone()).unwrap();
514    sigma.set(&[0, 0, 3, 2], cn1.clone()).unwrap();
515    sigma.set(&[0, 1, 0, 2], cn1.clone()).unwrap();
516    sigma.set(&[0, 2, 1, 3], cn1.clone()).unwrap();
517    sigma.set(&[1, 2, 0, 3], cn1.clone()).unwrap();
518    sigma.set(&[1, 2, 1, 1], cn1.clone()).unwrap();
519    sigma.set(&[1, 2, 2, 0], cn1.clone()).unwrap();
520    sigma.set(&[2, 0, 1, 2], cn1.clone()).unwrap();
521    sigma.set(&[2, 1, 2, 2], cn1.clone()).unwrap();
522    sigma.set(&[2, 2, 3, 3], cn1.clone()).unwrap();
523    sigma.set(&[3, 2, 0, 0], cn1.clone()).unwrap();
524    sigma.set(&[3, 2, 2, 3], cn1.clone()).unwrap();
525    sigma.set(&[3, 2, 3, 1], cn1.clone()).unwrap();
526    sigma.set(&[0, 0, 2, 3], cni.clone()).unwrap();
527    sigma.set(&[0, 0, 3, 1], cni.clone()).unwrap();
528    sigma.set(&[0, 1, 1, 3], cni.clone()).unwrap();
529    sigma.set(&[1, 0, 2, 1], cni.clone()).unwrap();
530    sigma.set(&[1, 3, 0, 1], cni.clone()).unwrap();
531    sigma.set(&[1, 3, 3, 0], cni.clone()).unwrap();
532    sigma.set(&[2, 0, 0, 3], cni.clone()).unwrap();
533    sigma.set(&[2, 0, 1, 1], cni.clone()).unwrap();
534    sigma.set(&[2, 1, 3, 3], cni.clone()).unwrap();
535    sigma.set(&[3, 0, 0, 1], cni.clone()).unwrap();
536    sigma.set(&[3, 3, 1, 0], cni.clone()).unwrap();
537    sigma.set(&[3, 3, 2, 1], cni.clone()).unwrap();
538
539    sigma
540}
541
542pub fn hep_lib<Aind: AbsInd, T: TensorLibraryData + Clone + Default>(
543    one: T,
544    zero: T,
545) -> TensorLibrary<MixedTensor<T, ExplicitKey<Aind>>, Aind>
546where
547{
548    let mut weyl = TensorLibrary::new();
549    initialize();
550    weyl.update_ids();
551
552    let gamma_key = PermutedStructure::identity(
553        gamma_data_weyl(AGS.gamma_strct::<Aind>(4), one.clone(), zero.clone()).into(),
554    );
555    // println!("permutation{}", gamma_key.rep_permutation);
556    weyl.insert_explicit(gamma_key);
557    let gamma_conj_key = PermutedStructure::identity(
558        gamma_conj_data_weyl(AGS.gamma_conj_strct::<Aind>(4), one.clone(), zero.clone()).into(),
559    );
560    // println!("permutation{}", gamma_key.rep_permutation);
561    weyl.insert_explicit(gamma_conj_key);
562    let gamma_adj_key = PermutedStructure::identity(
563        gamma_adj_data_weyl(AGS.gamma_adj_strct::<Aind>(4), one.clone(), zero.clone()).into(),
564    );
565    // println!("permutation{}", gamma_key.rep_permutation);
566    weyl.insert_explicit(gamma_adj_key);
567    let gamma0_key = PermutedStructure::identity(
568        gamma0_weyl(AGS.gamma0_strct::<Aind>(4), one.clone(), zero.clone()).into(),
569    );
570    // println!("permutation{}", gamma_key.rep_permutation);
571    weyl.insert_explicit(gamma0_key);
572
573    let gamma5_key = PermutedStructure::identity(
574        gamma5_weyl_data(AGS.gamma5_strct::<Aind>(4), one.clone(), zero.clone()).into(),
575    );
576    weyl.insert_explicit(gamma5_key);
577
578    let projm_key = PermutedStructure::identity(
579        proj_m_data_weyl(AGS.projm_strct::<Aind>(4), one.clone(), zero.clone()).into(),
580    );
581    weyl.insert_explicit(projm_key);
582
583    let projp_key = PermutedStructure::identity(
584        proj_p_data_weyl(AGS.projp_strct::<Aind>(4), one.clone(), zero.clone()).into(),
585    );
586    weyl.insert_explicit(projp_key);
587
588    weyl
589}
590
591pub fn insert_su3_color_tensors<Aind: AbsInd>(
592    lib: &mut TensorLibrary<MixedTensor<f64, ExplicitKey<Aind>>, Aind>,
593) {
594    initialize();
595
596    let t_key = PermutedStructure::identity(su3_generator_data(CS.t_strct::<Aind>(3, 8)).into());
597    lib.insert_explicit(t_key);
598
599    let f_key = PermutedStructure::identity(su3_structure_f_data(CS.f_strct::<Aind>(8)).into());
600    lib.insert_explicit(f_key);
601}
602
603pub fn hep_lib_su3<Aind: AbsInd>() -> TensorLibrary<MixedTensor<f64, ExplicitKey<Aind>>, Aind> {
604    let mut lib = hep_lib(1., 0.);
605    insert_su3_color_tensors(&mut lib);
606    lib
607}
608
609pub fn hep_lib_atom<Aind: AbsInd, T: TensorLibraryData + Clone + Default>()
610-> TensorLibrary<MixedTensor<T, ExplicitKey<Aind>>, Aind>
611where
612{
613    let mut weyl = TensorLibrary::new();
614    initialize();
615    weyl.update_ids();
616
617    let one = Atom::one();
618    let zero = Atom::Zero;
619
620    let gamma_key = PermutedStructure::identity(ParamOrConcrete::param(
621        gamma_data_weyl(AGS.gamma_strct::<Aind>(4), one.clone(), zero.clone())
622            .map_data(|a| a.re + a.im * Atom::i())
623            .into(),
624    ));
625    // println!("permutation{}", gamma_key.rep_permutation);
626    weyl.insert_explicit(gamma_key);
627    let gamma_conj_key = PermutedStructure::identity(ParamOrConcrete::param(
628        gamma_conj_data_weyl(AGS.gamma_conj_strct::<Aind>(4), one.clone(), zero.clone())
629            .map_data(|a| a.re + a.im * Atom::i())
630            .into(),
631    ));
632    // println!("permutation{}", gamma_key.rep_permutation);
633    weyl.insert_explicit(gamma_conj_key);
634    let gamma_adj_key = PermutedStructure::identity(ParamOrConcrete::param(
635        gamma_adj_data_weyl(AGS.gamma_adj_strct::<Aind>(4), one.clone(), zero.clone())
636            .map_data(|a| a.re + a.im * Atom::i())
637            .into(),
638    ));
639    // println!("permutation{}", gamma_key.rep_permutation);
640    weyl.insert_explicit(gamma_adj_key);
641    let gamma0_key = PermutedStructure::identity(ParamOrConcrete::param(
642        gamma0_weyl(AGS.gamma0_strct::<Aind>(4), one.clone(), zero.clone())
643            .map_data(|a| a.re + a.im * Atom::i())
644            .into(),
645    ));
646    // println!("permutation{}", gamma_key.rep_permutation);
647    weyl.insert_explicit(gamma0_key);
648
649    let gamma5_key = PermutedStructure::identity(ParamOrConcrete::param(
650        gamma5_weyl_data(AGS.gamma5_strct::<Aind>(4), one.clone(), zero.clone())
651            .map_data(|a| a.re + a.im * Atom::i())
652            .into(),
653    ));
654    weyl.insert_explicit(gamma5_key);
655
656    let projm_key = PermutedStructure::identity(ParamOrConcrete::param(
657        proj_m_data_weyl(AGS.projm_strct::<Aind>(4), one.clone(), zero.clone())
658            .map_data(|a| a.re + a.im * Atom::i())
659            .into(),
660    ));
661    weyl.insert_explicit(projm_key);
662
663    let projp_key = PermutedStructure::identity(ParamOrConcrete::param(
664        proj_p_data_weyl(AGS.projp_strct::<Aind>(4), one.clone(), zero.clone())
665            .map_data(|a| a.re + a.im * Atom::i())
666            .into(),
667    ));
668    weyl.insert_explicit(projp_key);
669
670    let color_t_key = PermutedStructure::identity(ParamOrConcrete::param(
671        su3_generator_data_atom(CS.t_strct::<Aind>(3, 8)).into(),
672    ));
673    weyl.insert_explicit(color_t_key);
674
675    let color_f_key = PermutedStructure::identity(ParamOrConcrete::param(
676        su3_structure_f_data_atom(CS.f_strct::<Aind>(8)).into(),
677    ));
678    weyl.insert_explicit(color_f_key);
679
680    weyl
681}
682
683pub type HepTensor<Aind> = MixedTensor<f64, ShadowedStructure<Aind>>;
684
685pub type HepNet<Aind> =
686    Network<NetworkStore<HepTensor<Aind>, Atom>, ExplicitKey<Aind>, Symbol, Aind>;
687
688pub static HEP_LIB: LazyLock<
689    TensorLibrary<MixedTensor<f64, ExplicitKey<AbstractIndex>>, AbstractIndex>,
690> = LazyLock::new(hep_lib_su3);
691
692pub static FUN_LIB: LazyLock<
693    SymbolLib<RealOrComplexTensor<f64, ShadowedStructure<AbstractIndex>>, PanicMissingConcrete>,
694> = LazyLock::new(|| {
695    let mut lib = PanicMissingConcrete::new_lib();
696    lib.insert(INBUILTS.conj, |a| match a {
697        RealOrComplexTensor::Complex(c) => RealOrComplexTensor::Complex(c.map_data(|x| x.conj())),
698        RealOrComplexTensor::Real(r) => RealOrComplexTensor::Real(r),
699    });
700    lib
701});
702
703#[cfg(test)]
704mod tests {
705
706    use spenso::{
707        network::{
708            Network, SingleSmallestDegree, SmallestDegreeIter, Steps,
709            parsing::{ParseSettings, ShadowedStructure, StrictTensorFilter},
710            store::NetworkStore,
711        },
712        structure::{HasStructure, abstract_index::AbstractIndex},
713    };
714    use symbolica::{
715        atom::{Atom, Symbol},
716        parse, parse_lit,
717    };
718
719    use super::*;
720
721    #[test]
722    fn simple_scalar() {
723        initialize();
724        let _a = HEP_LIB.get(&AGS.gamma_strct(4)).unwrap();
725
726        let expr = parse!("gamma(bis(4,l_5),bis(4,l_4),mink(4,l_4))*gamma(bis(4,l_6),bis(4,l_5),mink(4,l_4))*gamma(bis(4,l_4),bis(4,l_6),mink(4,l_5))*p(mink(4,l_5))
727            ",default_namespace="spenso");
728        // let expr = parse!(
729        // "gamma(bis(4,l_4),bis(4,l_6),mink(4,l_5))*p(mink(4,l_5))
730        // ",
731        // "spenso"
732        // );
733        // println!("{}", expr);
734
735        let mut net = Network::<
736            NetworkStore<MixedTensor<f64, ShadowedStructure<AbstractIndex>>, Atom>,
737            _,
738            Symbol,
739        >::try_from_view(
740            expr.as_view(),
741            &*HEP_LIB,
742            &ParseSettings::default().with_strict_tensor_filter(StrictTensorFilter::ContainsReps),
743        )
744        .unwrap();
745
746        println!(
747            "{}",
748            net.dot_display_impl(
749                |a| a.to_string(),
750                |a| Some(format!("{}", a.global_name.unwrap())),
751                |a| a.structure().global_name.unwrap().to_string(),
752                |a| a.to_string()
753            )
754        );
755
756        net.execute::<Steps<1>, SmallestDegreeIter<1>, _, _, _>(&*HEP_LIB, &*FUN_LIB)
757            .unwrap();
758        println!(
759            "{}",
760            net.dot_display_impl(
761                |a| a.to_string(),
762                |a| Some(format!("{}", a.global_name?)),
763                |a| a
764                    .structure()
765                    .global_name
766                    .map(|a| a.to_string())
767                    .unwrap_or("".to_string()),
768                |a| a.to_string()
769            )
770        );
771        net.execute::<Steps<1>, SmallestDegreeIter<2>, _, _, _>(&*HEP_LIB, &*FUN_LIB)
772            .unwrap();
773        println!(
774            "{}",
775            net.dot_display_impl(
776                |a| a.to_string(),
777                |a| Some(format!("{}", a.global_name?)),
778                |a| a
779                    .structure()
780                    .global_name
781                    .map(|a| a.to_string())
782                    .unwrap_or("".to_string()),
783                |a| a.to_string()
784            )
785        );
786
787        println!(
788            "{}",
789            net.dot_display_impl(
790                |a| a.to_string(),
791                |_| None,
792                |a| a.to_string(),
793                |a| a.to_string()
794            )
795        );
796        // if let ExecutionResult::Val(TensorOrScalarOrKey::Tensor { tensor, .. }) =
797        //     net.result().unwrap()
798        // {
799        //     // println!("YaY:{}", (&expr - &tensor.expression).expand());
800        //     // assert_eq!(expr, tensor.expression);
801        // } else {
802        //     panic!("Not tensor")
803        // }
804    }
805
806    #[test]
807    // #[should_panic]
808    fn parse_problem() {
809        initialize();
810        let _a = HEP_LIB.get(&AGS.gamma_strct(4)).unwrap();
811
812        let expr = parse_lit!(
813            (-1 * G
814                ^ 3 * P(0, mink(4, 0))
815                    * P(2, mink(4, 26))
816                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
817                    * gamma(bis(4, 7), bis(4, 6), mink(4, 1))
818                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
819                    + -1 * G
820                ^ 3 * P(0, mink(4, 26))
821                    * P(1, mink(4, 1))
822                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
823                    * gamma(bis(4, 7), bis(4, 6), mink(4, 0))
824                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
825                    + -1 * G
826                ^ 3 * P(0, mink(4, 26))
827                    * P(1, mink(4, 5))
828                    * g(mink(4, 0), mink(4, 1))
829                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
830                    * gamma(bis(4, 7), bis(4, 6), mink(4, 5))
831                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
832                    + -1 * G
833                ^ 3 * P(0, mink(4, 5))
834                    * P(2, mink(4, 26))
835                    * g(mink(4, 0), mink(4, 1))
836                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
837                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
838                    * gamma(bis(4, 7), bis(4, 6), mink(4, 5))
839                    + -1 * G
840                ^ 3 * P(1, mink(4, 1))
841                    * P(1, mink(4, 26))
842                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
843                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
844                    * gamma(bis(4, 7), bis(4, 6), mink(4, 0))
845                    + -1 * G
846                ^ 3 * P(1, mink(4, 26))
847                    * P(1, mink(4, 5))
848                    * g(mink(4, 0), mink(4, 1))
849                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
850                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
851                    * gamma(bis(4, 7), bis(4, 6), mink(4, 5))
852                    + -2 * G
853                ^ 3 * P(0, mink(4, 1))
854                    * P(0, mink(4, 26))
855                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
856                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
857                    * gamma(bis(4, 7), bis(4, 6), mink(4, 0))
858                    + -2 * G
859                ^ 3 * P(0, mink(4, 1))
860                    * P(1, mink(4, 26))
861                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
862                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
863                    * gamma(bis(4, 7), bis(4, 6), mink(4, 0))
864                    + -2 * G
865                ^ 3 * P(0, mink(4, 5))
866                    * Q(0, mink(4, 5))
867                    * g(mink(4, 0), mink(4, 1))
868                    * gamma(bis(4, 3), bis(4, 2), mink(4, 4))
869                    + -2 * G
870                ^ 3 * P(1, mink(4, 0))
871                    * P(1, mink(4, 1))
872                    * gamma(bis(4, 3), bis(4, 2), mink(4, 4))
873                    + -2 * G
874                ^ 3 * P(1, mink(4, 0))
875                    * P(2, mink(4, 26))
876                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
877                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
878                    * gamma(bis(4, 7), bis(4, 6), mink(4, 1))
879                    + -2 * G
880                ^ 3 * P(1, mink(4, 1))
881                    * P(2, mink(4, 0))
882                    * gamma(bis(4, 3), bis(4, 2), mink(4, 4))
883                    + -2 * G
884                ^ 3 * P(1, mink(4, 5))
885                    * P(2, mink(4, 5))
886                    * g(mink(4, 0), mink(4, 1))
887                    * gamma(bis(4, 3), bis(4, 2), mink(4, 4))
888                    + -4 * G
889                ^ 3 * P(0, mink(4, 1))
890                    * P(2, mink(4, 0))
891                    * gamma(bis(4, 3), bis(4, 2), mink(4, 4))
892                    + 2 * G
893                ^ 3 * P(0, mink(4, 0))
894                    * P(0, mink(4, 1))
895                    * gamma(bis(4, 3), bis(4, 2), mink(4, 4))
896                    + 2 * G
897                ^ 3 * P(0, mink(4, 0))
898                    * P(2, mink(4, 1))
899                    * gamma(bis(4, 3), bis(4, 2), mink(4, 4))
900                    + 2 * G
901                ^ 3 * P(0, mink(4, 1))
902                    * P(2, mink(4, 26))
903                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
904                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
905                    * gamma(bis(4, 7), bis(4, 6), mink(4, 0))
906                    + 2 * G
907                ^ 3 * P(0, mink(4, 26))
908                    * P(1, mink(4, 0))
909                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
910                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
911                    * gamma(bis(4, 7), bis(4, 6), mink(4, 1))
912                    + 2 * G
913                ^ 3 * P(0, mink(4, 5))
914                    * P(2, mink(4, 5))
915                    * g(mink(4, 0), mink(4, 1))
916                    * gamma(bis(4, 3), bis(4, 2), mink(4, 4))
917                    + 2 * G
918                ^ 3 * P(1, mink(4, 0))
919                    * P(1, mink(4, 26))
920                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
921                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
922                    * gamma(bis(4, 7), bis(4, 6), mink(4, 1))
923                    + 2 * G
924                ^ 3 * P(1, mink(4, i))
925                ^ 2 * g(mink(4, 0), mink(4, 1)) * gamma(bis(4, 3), bis(4, 2), mink(4, 4)) + G
926                ^ 3 * P(1, mink(4, 1))
927                    * P(2, mink(4, 26))
928                    * gamma(bis(4, 3), bis(4, 7), mink(4, 4))
929                    * gamma(bis(4, 6), bis(4, 2), mink(4, 26))
930                    * gamma(bis(4, 7), bis(4, 6), mink(4, 0))),
931            default_namespace = "spenso"
932        );
933        // println!("{}", expr);
934
935        let mut net = Network::<
936            NetworkStore<MixedTensor<f64, ShadowedStructure<AbstractIndex>>, Atom>,
937            _,
938            Symbol,
939        >::try_from_view(
940            expr.as_view(),
941            &*HEP_LIB,
942            &ParseSettings::default().with_strict_tensor_filter(StrictTensorFilter::ContainsReps),
943        )
944        .unwrap();
945
946        net.merge_ops();
947        println!(
948            "{}",
949            net.dot_display_impl(
950                |a| a.to_string(),
951                |a| Some(format!("{}", a.global_name.unwrap())),
952                |a| a.structure().global_name.unwrap().to_string(),
953                |a| a.to_string()
954            )
955        );
956
957        // net.validate();
958        net.execute::<Steps<1>, SingleSmallestDegree<true>, _, _, _>(&*HEP_LIB, &(*FUN_LIB))
959            .unwrap();
960        // net.validate();
961        // net.execute::<Steps<1>, SmallestDegree, _, _>(&*HEP_LIB);
962        // net.validate();
963        // net.execute::<Steps<1>, SmallestDegree, _, _>(&*HEP_LIB);
964        // net.validate();
965        // net.execute::<Steps<1>, SmallestDegree, _, _>(&*HEP_LIB);
966        // net.validate();
967        // net.execute::<Steps<1>, SmallestDegree, _, _>(&*HEP_LIB);
968        // net.validate();
969        // net.execute::<Steps<1>, SmallestDegree, _, _>(&*HEP_LIB);
970        // net.validate();
971        // net.execute::<Steps<1>, SmallestDegree, _, _>(&*HEP_LIB);
972        // net.validate();
973        // net.execute::<Steps<1>, SmallestDegree, _, _>(&*HEP_LIB);
974        // net.validate();
975        // net.execute::<Steps<1>, SmallestDegree, _, _>(&*HEP_LIB);
976        // net.validate();
977        // net.execute::<Steps<1>, SmallestDegree, _, _>(&*HEP_LIB);
978        // net.validate();
979        // net.execute::<Steps<1>, ContractScalars, _, _>(&*HEP_LIB);
980        // net.execute::<Steps<1>, SmallestDegree, _, _>(&*HEP_LIB);
981        // net.execute::<StepsDebug<1>, SingleSmallestDegree<true>, _, _>(&*HEP_LIB);
982        // net.execute::<Steps<1>, ContractScalars, _, _>(&*HEP_LIB);
983
984        //     .unwrap();
985        // net.execute::<Steps<14>, SingleSmallestDegree<false>, _, _>(&*HEP_LIB)
986        //     .unwrap();
987        // net.execute::<Steps<1>, SingleSmallestDegree<true>, _, _>(&*HEP_LIB)
988        //     .unwrap();
989        // // net.execute::<Sequential, SmallestDegree, _, _>(&*HEP_LIB)
990        //     .unwrap();
991        // println!(
992        //     "{}",
993        //     net.dot_display_impl(|a| a.to_string(), |_| None, |a| a.to_string())
994        // );
995
996        println!(
997            "{}",
998            net.dot_display_impl(
999                |a| a.to_string(),
1000                |_| None,
1001                |a| a.structure().to_string().replace('\n', "\\n"),
1002                |a| a.to_string()
1003            )
1004        );
1005        // if let ExecutionResult::Val(TensorOrScalarOrKey::Tensor { tensor, .. }) =
1006        //     net.result().unwrap()
1007        // {
1008        //     // println!("YaY:{}", (&expr - &tensor.expression).expand());
1009        //     // assert_eq!(expr, tensor.expression);
1010        // } else {
1011        //     panic!("Not tensor")
1012        // }
1013    }
1014
1015    #[test]
1016    fn transpose_test() {}
1017}
spenso_hep_lib/lib.rs

spenso_hep_lib/
lib.rs
