momtrop/lib.rs
1//! `momtrop` is a flexible Rust library for tropical Feynman sampling of loop integrals in momentum space (https://arxiv.org/abs/2504.09613).
2//! It allows users to define graphs with arbitrary edge weights, dimensions, and external kinematics,
3//! and supports custom floating-point types via the `MomtropFloat` trait.
4//!
5//! The library only takes care of sample generation, and defers the evaluation of the integrand to the user.
6
7use std::iter::repeat_with;
8
9use bincode::{Decode, Encode};
10
11use float::MomTropFloat;
12use itertools::Itertools;
13use matrix::{DecompositionResult, SquareMatrix};
14use preprocessing::{TropicalGraph, TropicalSubgraphTable};
15use rand::Rng;
16use sampling::{SamplingError, sample};
17use serde::{Deserialize, Serialize};
18use vector::Vector;
19
20pub mod float;
21pub mod gamma;
22pub mod matrix;
23mod mimic_rng;
24mod preprocessing;
25mod sampling;
26pub mod vector;
27
28#[cfg(feature = "python_api")]
29pub mod python;
30
31/// Maximum number of edges supported by momtrop.
32pub const MAX_EDGES: usize = 64;
33/// Maximum number of vertices supported by momtrop.
34pub const MAX_VERTICES: usize = 256;
35
36#[derive(Debug, Clone)]
37/// Struct containing all runtime settings.
38pub struct TropicalSamplingSettings {
39 /// `matrix_stability_test` tests the numerical stability
40 /// of some of the matrix routines used during sampling. This is done by checking
41 /// how far L L^-1 is from the identity matrix in terms of a L_2_1 norm. If this distance is
42 /// larger than the float provided, the sampling will return an Err.
43 pub matrix_stability_test: Option<f64>,
44 /// Setting `print_debug_info` to `true` will print out the results of intermediate steps during sampling.
45 /// Useful for debugging.
46 pub print_debug_info: bool,
47 /// Enabling `return_metadata` provides access to some of the intermediate results. This is useful for
48 /// exploring new applications.
49 pub return_metadata: bool,
50}
51
52#[allow(clippy::derivable_impls)]
53impl Default for TropicalSamplingSettings {
54 fn default() -> Self {
55 Self {
56 matrix_stability_test: None,
57 print_debug_info: false,
58 return_metadata: false,
59 }
60 }
61}
62
63pub fn assert_approx_eq<T: MomTropFloat>(res: &T, target: &T, tolerance: &T) {
64 if approx_eq(res, target, tolerance) {
65 } else {
66 panic!(
67 "assert_approx_eq failed: \n{:?} != \n{:?} with tolerance {:?}",
68 &res, &target, &tolerance
69 )
70 }
71}
72
73/// Main graph struct from which a sampler can be build. This graph should be stripped of
74/// tree-like and external edges.
75/// Vertices which would have external edges attached to them must be added to the `externals` field.
76#[derive(Clone)]
77pub struct Graph {
78 pub edges: Vec<Edge>,
79 pub externals: Vec<u8>,
80}
81
82/// Edge struct from which a graph can be specified.
83#[derive(Clone, Copy)]
84pub struct Edge {
85 pub vertices: (u8, u8),
86 pub is_massive: bool,
87 /// `weight` specifies the power of the propagator corresponding to an edge.
88 pub weight: f64,
89}
90
91#[derive(Debug, Clone)]
92/// Return type of the sample function
93pub struct TropicalSampleResult<T: MomTropFloat, const D: usize> {
94 /// loop momenta of sample point.
95 pub loop_momenta: Vec<Vector<T, D>>,
96 /// tropical approximation of U at sample point.
97 pub u_trop: T,
98 /// tropical approximation of V at sample point.
99 pub v_trop: T,
100 /// U at sample point.
101 pub u: T,
102 /// V at sample point.
103 pub v: T,
104 /// Jacobian at sample point.
105 pub jacobian: T,
106 /// metadata, is `None` is `return_metadata` is disabled.
107 pub metadata: Option<Metadata<T, D>>,
108}
109
110#[derive(Debug, Clone)]
111/// Optional metadata for advanced users, see \citation for the definition of the fields
112/// in this struct.
113pub struct Metadata<T: MomTropFloat, const D: usize> {
114 /// Vectors q sampled from Gaussian.
115 pub q_vectors: Vec<Vector<T, D>>,
116 /// Parameter lambda sampled from Gamma distribution.
117 pub lambda: T,
118 /// Matrix L at sample point.
119 pub l_matrix: SquareMatrix<T>,
120 /// Cholesky decomposition and inverse of L at sample point.
121 pub decompoisiton_result: DecompositionResult<T>,
122 /// The vectors u at sample point.
123 pub u_vectors: Vec<Vector<T, D>>,
124 /// L^-1 u at sample point
125 pub shift: Vec<Vector<T, D>>,
126}
127
128impl Graph {
129 /// Build the sampler from a graph. `loop_signature` is the signature matrix determining the loop-momentum routing.
130 /// Returns `Err` if a divergent subgraph is found.
131 ///
132 /// For example, the following code builds a sampler for the massless triangle with all edges having weight 2/3.
133 ///
134 ///```rust
135 ///use momtrop::{Graph, Edge};
136 ///
137 ///let weight = 2.0 / 3.0;
138 ///
139 ///let triangle_edges = vec![
140 /// Edge {
141 /// vertices: (0, 1),
142 /// is_massive: false,
143 /// weight,
144 /// },
145 /// Edge {
146 /// vertices: (1, 2),
147 /// is_massive: false,
148 /// weight,
149 /// },
150 /// Edge {
151 /// vertices: (2, 0),
152 /// is_massive: false,
153 /// weight,
154 /// },
155 ///];
156 ///
157 ///let externals = vec![0, 1, 2];
158 ///
159 ///let triangle_graph = Graph {
160 /// edges: triangle_edges,
161 /// externals,
162 ///};
163 ///
164 /// let loop_signature = vec![vec![1]; 3];
165 /// let triangle_sampler = triangle_graph.build_sampler::<3>(loop_signature);
166 ///
167 /// assert!(triangle_sampler.is_ok());
168 /// ```
169 pub fn build_sampler<const D: usize>(
170 self,
171 loop_signature: Vec<Vec<isize>>,
172 ) -> Result<SampleGenerator<D>, String> {
173 let tropical_graph = TropicalGraph::from_graph(self, D)?;
174 let table = TropicalSubgraphTable::generate_from_tropical(tropical_graph, D, None)?;
175
176 Ok(SampleGenerator {
177 loop_signature,
178 table,
179 })
180 }
181}
182
183fn approx_eq<T: MomTropFloat>(res: &T, target: &T, tolerance: &T) -> bool {
184 if target == &res.zero() {
185 &res.abs() < tolerance
186 } else {
187 &((res.ref_sub(target)) / target).abs() < tolerance
188 }
189}
190
191#[derive(Clone, Debug, Serialize, Deserialize, Decode, Encode)]
192/// Sampler struct from which sample points can be generated.
193pub struct SampleGenerator<const D: usize> {
194 loop_signature: Vec<Vec<isize>>,
195 table: TropicalSubgraphTable,
196}
197
198impl<const D: usize> SampleGenerator<D> {
199 /// Generate a sample point for a given random point `x_space_point` in the unit hypercube.
200 /// `edge_data` contains masses and shifs of each propagator.
201 /// Dimensionality of the unit hypercube, should match the length of `x_space_point`.
202 ///```rust
203 ///# use momtrop::{Graph, Edge};
204 ///use momtrop::{TropicalSamplingSettings, vector::Vector};
205 ///use rand::{rngs::StdRng, SeedableRng, Rng};
206 ///use std::iter::repeat_with;
207 ///#
208 ///# let weight = 2.0 / 3.0;
209 ///#
210 ///# let triangle_edges = vec![
211 ///# Edge {
212 ///# vertices: (0, 1),
213 ///# is_massive: false,
214 ///# weight,
215 ///# },
216 ///# Edge {
217 ///# vertices: (1, 2),
218 ///# is_massive: false,
219 ///# weight,
220 ///# },
221 ///# Edge {
222 ///# vertices: (2, 0),
223 ///# is_massive: false,
224 ///# weight,
225 ///# },
226 ///# ];
227 ///#
228 ///# let externals = vec![0, 1, 2];
229 ///#
230 ///# let triangle_graph = Graph {
231 ///# edges: triangle_edges,
232 ///# externals,
233 ///# };
234 ///#
235 ///# let loop_signature = vec![vec![1]; 3];
236 ///let triangle_sampler = triangle_graph.build_sampler::<3>(loop_signature).unwrap();
237 ///
238 ///let settings = TropicalSamplingSettings::default();
239 ///
240 ///let p_1 = Vector::from_array([0.1, 0.2, 0.3]);
241 ///let p_2 = Vector::from_array([0.4, 0.5, 0.6]);
242 ///let edge_data = vec![(None, Vector::new_from_num(&0.0)), (None, p_1), (None, p_2)];
243 ///
244 ///let mut rng: StdRng = SeedableRng::seed_from_u64(42);
245 ///
246 ///let x_space_point = repeat_with(|| rng.r#gen::<f64>()).take(triangle_sampler.get_dimension()).collect::<Vec<_>>();
247 ///let sample = triangle_sampler.generate_sample_from_x_space_point(&x_space_point, edge_data, &settings, None);
248 /// assert!(sample.is_ok());
249 ///
250 /// ```
251 pub fn generate_sample_from_x_space_point<T: MomTropFloat>(
252 &self,
253 x_space_point: &[T],
254 edge_data: Vec<(Option<T>, vector::Vector<T, D>)>,
255 settings: &TropicalSamplingSettings,
256 force_sector: Option<&[usize]>,
257 ) -> Result<TropicalSampleResult<T, D>, SamplingError> {
258 sample(
259 &self.table,
260 x_space_point,
261 &self.loop_signature,
262 &edge_data,
263 settings,
264 force_sector,
265 )
266 }
267
268 /// Alternative to `generate_sample_from_x_space_point`. Uses a `Rng` to generate the point
269 /// in the unit hypercube.
270 ///```rust
271 ///# use momtrop::{Graph, Edge};
272 ///use momtrop::{TropicalSamplingSettings, vector::Vector};
273 ///use rand::{rngs::StdRng, SeedableRng, Rng};
274 ///use std::iter::repeat_with;
275 ///#
276 ///# let weight = 2.0 / 3.0;
277 ///#
278 ///# let triangle_edges = vec![
279 ///# Edge {
280 ///# vertices: (0, 1),
281 ///# is_massive: false,
282 ///# weight,
283 ///# },
284 ///# Edge {
285 ///# vertices: (1, 2),
286 ///# is_massive: false,
287 ///# weight,
288 ///# },
289 ///# Edge {
290 ///# vertices: (2, 0),
291 ///# is_massive: false,
292 ///# weight,
293 ///# },
294 ///# ];
295 ///#
296 ///# let externals = vec![0, 1, 2];
297 ///#
298 ///# let triangle_graph = Graph {
299 ///# edges: triangle_edges,
300 ///# externals,
301 ///# };
302 ///#
303 ///# let loop_signature = vec![vec![1]; 3];
304 ///let triangle_sampler = triangle_graph.build_sampler::<3>(loop_signature).unwrap();
305 ///
306 ///let settings = TropicalSamplingSettings::default();
307 ///
308 ///let p_1 = Vector::from_array([0.1, 0.2, 0.3]);
309 ///let p_2 = Vector::from_array([0.4, 0.5, 0.6]);
310 ///let edge_data = vec![(None, Vector::new_from_num(&0.0)), (None, p_1), (None, p_2)];
311 ///
312 ///let mut rng: StdRng = SeedableRng::seed_from_u64(42);
313 ///
314 ///let sample = triangle_sampler.generate_sample_from_rng(edge_data, &settings, &mut rng);
315 /// assert!(sample.is_ok());
316 ///
317 /// ```
318 pub fn generate_sample_from_rng<T: MomTropFloat, R: Rng>(
319 &self,
320 edge_data: Vec<(Option<T>, vector::Vector<T, D>)>,
321 settings: &TropicalSamplingSettings,
322 rng: &mut R,
323 ) -> Result<TropicalSampleResult<T, D>, SamplingError> {
324 let const_builder = edge_data[0].1.zero();
325
326 let num_vars = self.get_dimension();
327 let x_space_point = repeat_with(|| const_builder.from_f64(rng.r#gen::<f64>()))
328 .take(num_vars)
329 .collect_vec();
330
331 self.generate_sample_from_x_space_point(&x_space_point, edge_data, settings, None)
332 }
333
334 /// Dimensionality of the unit hypercube, should match the length of `x_space_point`.
335 ///```rust
336 ///# use momtrop::{Graph, Edge};
337 ///#
338 ///# let weight = 2.0 / 3.0;
339 ///#
340 ///# let triangle_edges = vec![
341 ///# Edge {
342 ///# vertices: (0, 1),
343 ///# is_massive: false,
344 ///# weight,
345 ///# },
346 ///# Edge {
347 ///# vertices: (1, 2),
348 ///# is_massive: false,
349 ///# weight,
350 ///# },
351 ///# Edge {
352 ///# vertices: (2, 0),
353 ///# is_massive: false,
354 ///# weight,
355 ///# },
356 ///# ];
357 ///#
358 ///# let externals = vec![0, 1, 2];
359 ///#
360 ///# let triangle_graph = Graph {
361 ///# edges: triangle_edges,
362 ///# externals,
363 ///# };
364 ///#
365 ///# let loop_signature = vec![vec![1]; 3];
366 ///let triangle_sampler = triangle_graph.build_sampler::<3>(loop_signature).unwrap();
367 ///
368 /// assert_eq!(triangle_sampler.get_dimension(), 2*3 - 1 + 3 + (3 % 2) );
369 /// ```
370 pub fn get_dimension(&self) -> usize {
371 self.table.get_num_variables()
372 }
373
374 /// Degree of divergence of the provided graph.
375 ///```rust
376 ///# use momtrop::{Graph, Edge};
377 ///#
378 ///# let weight = 2.0 / 3.0;
379 ///#
380 ///# let triangle_edges = vec![
381 ///# Edge {
382 ///# vertices: (0, 1),
383 ///# is_massive: false,
384 ///# weight,
385 ///# },
386 ///# Edge {
387 ///# vertices: (1, 2),
388 ///# is_massive: false,
389 ///# weight,
390 ///# },
391 ///# Edge {
392 ///# vertices: (2, 0),
393 ///# is_massive: false,
394 ///# weight,
395 ///# },
396 ///# ];
397 ///#
398 ///# let externals = vec![0, 1, 2];
399 ///#
400 ///# let triangle_graph = Graph {
401 ///# edges: triangle_edges,
402 ///# externals,
403 ///# };
404 ///#
405 ///# let loop_signature = vec![vec![1]; 3];
406 ///let triangle_sampler = triangle_graph.build_sampler::<3>(loop_signature).unwrap();
407 ///
408 /// assert_eq!(triangle_sampler.get_dod(), 3.0 * 2.0 / 3.0 - 3.0 / 2.0);
409 /// ```
410 pub fn get_dod(&self) -> f64 {
411 self.table.tropical_graph.dod
412 }
413
414 /// Iterate over the edge weights of the provided graph.
415 ///```rust
416 ///# use momtrop::{Graph, Edge};
417 ///#
418 ///# let weight = 2.0 / 3.0;
419 ///#
420 ///# let triangle_edges = vec![
421 ///# Edge {
422 ///# vertices: (0, 1),
423 ///# is_massive: false,
424 ///# weight,
425 ///# },
426 ///# Edge {
427 ///# vertices: (1, 2),
428 ///# is_massive: false,
429 ///# weight,
430 ///# },
431 ///# Edge {
432 ///# vertices: (2, 0),
433 ///# is_massive: false,
434 ///# weight,
435 ///# },
436 ///# ];
437 ///#
438 ///# let externals = vec![0, 1, 2];
439 ///#
440 ///# let triangle_graph = Graph {
441 ///# edges: triangle_edges,
442 ///# externals,
443 ///# };
444 ///#
445 ///# let loop_signature = vec![vec![1]; 3];
446 ///let triangle_sampler = triangle_graph.build_sampler::<3>(loop_signature).unwrap();
447 ///
448 /// assert!(triangle_sampler.iter_edge_weights().all(|w| w == 2.0/3.0));
449 /// ```
450 pub fn iter_edge_weights(&self) -> impl Iterator<Item = f64> + '_ {
451 self.table.tropical_graph.topology.iter().map(|e| e.weight)
452 }
453
454 /// Get the number of edges of the provided graph.
455 ///```rust
456 ///# use momtrop::{Graph, Edge};
457 ///#
458 ///# let weight = 2.0 / 3.0;
459 ///#
460 ///# let triangle_edges = vec![
461 ///# Edge {
462 ///# vertices: (0, 1),
463 ///# is_massive: false,
464 ///# weight,
465 ///# },
466 ///# Edge {
467 ///# vertices: (1, 2),
468 ///# is_massive: false,
469 ///# weight,
470 ///# },
471 ///# Edge {
472 ///# vertices: (2, 0),
473 ///# is_massive: false,
474 ///# weight,
475 ///# },
476 ///# ];
477 ///#
478 ///# let externals = vec![0, 1, 2];
479 ///#
480 ///# let triangle_graph = Graph {
481 ///# edges: triangle_edges,
482 ///# externals,
483 ///# };
484 ///#
485 ///# let loop_signature = vec![vec![1]; 3];
486 ///let triangle_sampler = triangle_graph.build_sampler::<3>(loop_signature).unwrap();
487 ///
488 /// assert_eq!(triangle_sampler.get_num_edges(), 3);
489 /// ```
490 pub fn get_num_edges(&self) -> usize {
491 self.table.tropical_graph.topology.len()
492 }
493
494 /// Get the smallest degree of divergence in the subgraph table.
495 pub fn get_smallest_dod(&self) -> f64 {
496 self.table.get_smallest_dod()
497 }
498}