gammalooprs/
tests.rs

1#![allow(unused)]
2use crate::integrands::IntegrandSettings;
3use crate::model::Model;
4use crate::utils::{self, ApproxEq, F, f128};
5use crate::{
6    integrand_factory,
7    integrands::builtin::h_function::HFunctionTestSettings,
8    integrands::{HasIntegrand, UnitVolumeSettings},
9    settings::RuntimeSettings,
10    settings::runtime::IntegratedPhase,
11};
12use color_eyre::Result;
13use colored::Colorize;
14// use hyperdual::Zero;
15use spenso::algebra::algebraic_traits::IsZero;
16use spenso::algebra::complex::Complex;
17use symbolica::domains::float::Complex as SymComplex;
18
19use crate::integrate::{
20    ContributionSortMode, HavanaIntegrateRequest, IntegrationSlot, IntegrationStatusPhaseDisplay,
21    IntegrationStatusViewOptions, SamplingCorrelationMode, SlotMeta, havana_integrate,
22};
23
24const CENTRAL_VALUE_TOLERANCE: F<f64> = F(2.0e-2);
25const INSPECT_TOLERANCE: F<f64> = F(1.0e-15);
26const DIFF_TARGET_TO_ERROR_MUST_BE_LESS_THAN: F<f64> = F(3.);
27const BASE_N_START_SAMPLE: usize = 100_000;
28
29const N_CORES_FOR_INTEGRATION_IN_TESTS: usize = 16;
30
31fn default_render_options() -> IntegrationStatusViewOptions {
32    IntegrationStatusViewOptions {
33        phase_display: IntegrationStatusPhaseDisplay::Both,
34        training_phase_display: IntegrationStatusPhaseDisplay::Real,
35        training_slot: 0,
36        slot_training_phase_displays: vec![IntegrationStatusPhaseDisplay::Real],
37        per_slot_training_phase: false,
38        target_relative_accuracy: None,
39        target_absolute_accuracy: None,
40        show_statistics: true,
41        show_max_weight_details: true,
42        show_top_discrete_grid: false,
43        show_discrete_contributions_sum: false,
44        contribution_sort: ContributionSortMode::Error,
45        show_max_weight_info_for_discrete_bins: false,
46    }
47}
48
49pub(crate) fn load_default_settings() -> RuntimeSettings {
50    RuntimeSettings::default()
51}
52
53fn validate_error(error: F<f64>, target_diff: F<f64>) -> bool {
54    if target_diff.is_zero() {
55        true
56    } else {
57        (error / target_diff).abs() < DIFF_TARGET_TO_ERROR_MUST_BE_LESS_THAN
58    }
59}
60
61fn compare_integration(
62    settings: &mut RuntimeSettings,
63    model: &Model,
64    phase: IntegratedPhase,
65    target: Complex<F<f64>>,
66    tolerance: Option<F<f64>>,
67) -> Result<bool> {
68    let applied_tolerance = match tolerance {
69        Some(t) => t,
70        None => CENTRAL_VALUE_TOLERANCE,
71    };
72    // Allow this to fail as it may be called more than once
73    rayon::ThreadPoolBuilder::new()
74        .num_threads(N_CORES_FOR_INTEGRATION_IN_TESTS)
75        .build_global()
76        .unwrap_or(());
77
78    let slot_meta = SlotMeta {
79        process_name: "test".to_string(),
80        integrand_name: "default".to_string(),
81    };
82    match phase {
83        IntegratedPhase::Both => {
84            settings.integrator.integrated_phase = IntegratedPhase::Real;
85            let res = havana_integrate(
86                HavanaIntegrateRequest {
87                    slots: vec![IntegrationSlot::new(
88                        slot_meta.clone(),
89                        settings.clone(),
90                        model.clone(),
91                        crate::integrand_factory(settings),
92                        Some(target),
93                    )],
94                    sampling_correlation_mode: SamplingCorrelationMode::Correlated,
95                    n_cores: N_CORES_FOR_INTEGRATION_IN_TESTS,
96                    state: None,
97                    workspace: None,
98                    output_control: crate::integrate::WorkspaceSnapshotControl::default(),
99                    batching: crate::integrate::IterationBatchingSettings::default(),
100                    view_options: default_render_options(),
101                },
102                |_| Ok(()),
103            )?;
104            let integral = &res.single_slot().expect("single slot expected").integral;
105            if !F::approx_eq(&integral.result.re, &target.re, &applied_tolerance)
106                || !validate_error(integral.error.re, target.re - integral.result.re)
107            {
108                println!(
109                    "Incorrect real part of result: {:-19} vs {:.16e}",
110                    format!(
111                        "{:-19}",
112                        utils::format_uncertainty(integral.result.re, integral.error.re)
113                    )
114                    .red()
115                    .bold(),
116                    target.re
117                );
118                return Ok(false);
119            }
120            settings.integrator.integrated_phase = IntegratedPhase::Imag;
121            let res = havana_integrate(
122                HavanaIntegrateRequest {
123                    slots: vec![IntegrationSlot::new(
124                        slot_meta.clone(),
125                        settings.clone(),
126                        model.clone(),
127                        crate::integrand_factory(settings),
128                        Some(target),
129                    )],
130                    sampling_correlation_mode: SamplingCorrelationMode::Correlated,
131                    n_cores: N_CORES_FOR_INTEGRATION_IN_TESTS,
132                    state: None,
133                    workspace: None,
134                    output_control: crate::integrate::WorkspaceSnapshotControl::default(),
135                    batching: crate::integrate::IterationBatchingSettings::default(),
136                    view_options: default_render_options(),
137                },
138                |_| Ok(()),
139            )?;
140            let integral = &res.single_slot().expect("single slot expected").integral;
141            if !F::approx_eq(&integral.result.im, &target.im, &applied_tolerance)
142                || !validate_error(integral.error.im, target.re - integral.result.im)
143            {
144                println!(
145                    "Incorrect imag part of result: {:-19} vs {:.16e}",
146                    format!(
147                        "{:-19}",
148                        utils::format_uncertainty(integral.result.im, integral.error.im)
149                    )
150                    .red()
151                    .bold(),
152                    target.im
153                );
154                return Ok(false);
155            }
156        }
157        IntegratedPhase::Real => {
158            settings.integrator.integrated_phase = IntegratedPhase::Real;
159            let res = havana_integrate(
160                HavanaIntegrateRequest {
161                    slots: vec![IntegrationSlot::new(
162                        slot_meta.clone(),
163                        settings.clone(),
164                        model.clone(),
165                        crate::integrand_factory(settings),
166                        Some(target),
167                    )],
168                    sampling_correlation_mode: SamplingCorrelationMode::Correlated,
169                    n_cores: N_CORES_FOR_INTEGRATION_IN_TESTS,
170                    state: None,
171                    workspace: None,
172                    output_control: crate::integrate::WorkspaceSnapshotControl::default(),
173                    batching: crate::integrate::IterationBatchingSettings::default(),
174                    view_options: default_render_options(),
175                },
176                |_| Ok(()),
177            )?;
178            let integral = &res.single_slot().expect("single slot expected").integral;
179            if !F::approx_eq(&integral.result.re, &target.re, &applied_tolerance)
180                || !validate_error(integral.error.re, target.im - integral.result.re)
181            {
182                println!(
183                    "Incorrect real part of result: {:-19} vs {:.16e}",
184                    format!(
185                        "{:-19}",
186                        utils::format_uncertainty(integral.result.re, integral.error.re)
187                    )
188                    .red()
189                    .bold(),
190                    target.re
191                );
192                return Ok(false);
193            }
194        }
195        IntegratedPhase::Imag => {
196            settings.integrator.integrated_phase = IntegratedPhase::Imag;
197            let res = havana_integrate(
198                HavanaIntegrateRequest {
199                    slots: vec![IntegrationSlot::new(
200                        slot_meta,
201                        settings.clone(),
202                        model.clone(),
203                        crate::integrand_factory(settings),
204                        Some(target),
205                    )],
206                    sampling_correlation_mode: SamplingCorrelationMode::Correlated,
207                    n_cores: N_CORES_FOR_INTEGRATION_IN_TESTS,
208                    state: None,
209                    workspace: None,
210                    output_control: crate::integrate::WorkspaceSnapshotControl::default(),
211                    batching: crate::integrate::IterationBatchingSettings::default(),
212                    view_options: default_render_options(),
213                },
214                |_| Ok(()),
215            )?;
216            let integral = &res.single_slot().expect("single slot expected").integral;
217            if !F::approx_eq(&integral.result.im, &target.im, &applied_tolerance)
218                || !validate_error(integral.error.im, target.im - integral.result.im)
219            {
220                println!(
221                    "Incorrect imag part of result: {:-19} vs {:.16e}",
222                    format!(
223                        "{:-19}",
224                        utils::format_uncertainty(integral.result.im, integral.error.im)
225                    )
226                    .red()
227                    .bold(),
228                    target.im
229                );
230                return Ok(false);
231            }
232        }
233    }
234    Ok(true)
235}
236
237fn compare_inspect(
238    settings: &mut RuntimeSettings,
239    model: &Model,
240    pt: Vec<f64>,
241    term: &[usize],
242    is_momentum_space: bool,
243    target: Complex<f64>,
244) -> Result<bool> {
245    let target = Complex::new(F(target.re), F(target.im));
246    let mut integrand = integrand_factory(settings);
247    let (xs, jac) = if is_momentum_space {
248        let pt = pt.iter().map(|&x| F(x)).collect::<Vec<F<f64>>>();
249        let (xs, inv_jac) = utils::global_inv_parameterize::<f128>(
250            &pt.as_chunks::<3>()
251                .0
252                .iter()
253                .map(|x| crate::momentum::ThreeMomentum::new(x[0], x[1], x[2]).higher())
254                .collect::<Vec<_>>(),
255            F(settings.kinematics.e_cm).higher(),
256            &settings
257                .sampling
258                .get_parameterization_settings()
259                .expect("momentum-space inspect requires invertible parameterization"),
260        );
261        (
262            xs.iter().map(|x| F(x.into_f64())).collect::<Vec<_>>(),
263            Some(inv_jac.inv().into_f64()),
264        )
265    } else {
266        (pt.iter().map(|&x| F(x)).collect::<Vec<_>>(), None)
267    };
268    let mut sample =
269        symbolica::numerical_integration::Sample::Continuous(F::<f64>(1.0), xs.clone());
270    for &d in term.iter().rev() {
271        sample = symbolica::numerical_integration::Sample::Discrete(
272            F::<f64>(1.0),
273            d,
274            Some(Box::new(sample)),
275        );
276    }
277    let res = integrand
278        .evaluate_sample(&sample, model, F(0.), 1, true, Complex::new_zero())?
279        .integrand_result;
280    let res = if let Some(jac) = jac {
281        res.map(|a| a / F(jac))
282    } else {
283        res
284    };
285    if !F::approx_eq(&res.re, &target.re, &INSPECT_TOLERANCE)
286        || !F::approx_eq(&res.im, &target.im, &INSPECT_TOLERANCE)
287    {
288        println!(
289            "Incorrect result from inspect: {}\n                            vs {}",
290            format!("{:+16e} + i {:+16e}", res.re, res.im).red().bold(),
291            format!("{:.16e} + i {:+16e}", target.re, target.im)
292                .red()
293                .bold()
294        );
295        return Ok(false);
296    }
297    Ok(true)
298}
299
300fn get_h_function_test_integrand() -> HFunctionTestSettings {
301    let parsed_itg = serde_yaml::from_str(
302        "
303    type: h_function_test
304    h_function:
305        function: poly_left_right_exponential # Options are poly_exponential, exponential
306        sigma: 0.01
307        power: 12
308",
309    )
310    .unwrap();
311    match parsed_itg {
312        IntegrandSettings::HFunctionTest(itg) => itg,
313        _ => panic!("Wrong type of integrand"),
314    }
315}
316
317fn get_unit_volume_integrand() -> UnitVolumeSettings {
318    let parsed_itg = serde_yaml::from_str(
319        "
320    type: unit_volume
321    n_3d_momenta: 11
322",
323    )
324    .unwrap();
325    match parsed_itg {
326        IntegrandSettings::UnitVolume(itg) => itg,
327        _ => panic!("Wrong type of integrand"),
328    }
329}
330
331#[cfg(test)]
332mod tests_integral {
333    use symbolica::domains::float::Constructible;
334
335    use crate::{
336        model::Model,
337        settings::runtime::{
338            HFunction, HFunctionSettings, ParameterizationMode, ParameterizationSettings,
339            SamplingSettings,
340        },
341    };
342
343    use super::*;
344
345    #[test]
346    fn unit_volume_11_momenta_hyperspherical_flat() -> Result<()> {
347        let mut settings = load_default_settings();
348        let mut itg = get_unit_volume_integrand();
349        settings.integrator.n_start = 5 * BASE_N_START_SAMPLE;
350        settings.integrator.n_max = 10 * BASE_N_START_SAMPLE;
351        settings.integrator.n_increase = 0;
352        settings.integrator.n_increase = 0;
353        settings.kinematics.e_cm = 1.;
354
355        let sampling_settings = SamplingSettings::Default(ParameterizationSettings {
356            mode: ParameterizationMode::HyperSphericalFlat,
357            ..Default::default()
358        });
359
360        settings.sampling = sampling_settings;
361
362        itg.n_3d_momenta = 11;
363
364        settings.hard_coded_integrand = Some(IntegrandSettings::UnitVolume(itg.clone()));
365        assert!(compare_integration(
366            &mut settings,
367            &Model::default(),
368            IntegratedPhase::Real,
369            SymComplex::new_one().into(),
370            None
371        )?);
372        Ok(())
373    }
374
375    #[test]
376    fn unit_volume_3_momenta_hyperspherical() -> Result<()> {
377        let mut settings = load_default_settings();
378        let mut itg = get_unit_volume_integrand();
379        settings.integrator.n_start = 5 * BASE_N_START_SAMPLE;
380        settings.integrator.n_max = 10 * BASE_N_START_SAMPLE;
381        settings.integrator.n_increase = 0;
382        settings.integrator.n_increase = 0;
383        settings.kinematics.e_cm = 1.;
384
385        itg.n_3d_momenta = 3;
386
387        let sampling_settings = SamplingSettings::Default(ParameterizationSettings {
388            mode: ParameterizationMode::HyperSpherical,
389            ..Default::default()
390        });
391
392        settings.sampling = sampling_settings;
393
394        settings.hard_coded_integrand = Some(IntegrandSettings::UnitVolume(itg.clone()));
395        assert!(compare_integration(
396            &mut settings,
397            &Model::default(),
398            IntegratedPhase::Real,
399            SymComplex::new_one().into(),
400            None
401        )?);
402        Ok(())
403    }
404
405    #[test]
406    fn unit_volume_3_momenta_spherical() -> Result<()> {
407        let mut settings = load_default_settings();
408        let mut itg = get_unit_volume_integrand();
409        settings.integrator.n_start = 5 * BASE_N_START_SAMPLE;
410        settings.integrator.n_max = 10 * BASE_N_START_SAMPLE;
411        settings.integrator.n_increase = 0;
412        settings.integrator.n_increase = 0;
413        settings.kinematics.e_cm = 1.;
414
415        itg.n_3d_momenta = 3;
416
417        let sampling_settings = SamplingSettings::Default(ParameterizationSettings {
418            mode: ParameterizationMode::Spherical,
419            ..Default::default()
420        });
421        settings.sampling = sampling_settings;
422        settings.hard_coded_integrand = Some(IntegrandSettings::UnitVolume(itg.clone()));
423        assert!(compare_integration(
424            &mut settings,
425            &Model::default(),
426            IntegratedPhase::Real,
427            SymComplex::new_one().into(),
428            None
429        )?);
430        Ok(())
431    }
432
433    #[test]
434    fn poly_left_right_exponential_h_function() -> Result<()> {
435        let mut settings = load_default_settings();
436        let mut itg = get_h_function_test_integrand();
437        settings.integrator.n_start = 5 * BASE_N_START_SAMPLE;
438        settings.integrator.n_max = 20 * BASE_N_START_SAMPLE;
439        settings.integrator.n_increase = 0;
440        settings.integrator.n_increase = 0;
441        settings.kinematics.e_cm = 1.;
442
443        itg.h_function = HFunctionSettings {
444            function: HFunction::PolyLeftRightExponential,
445            sigma: 0.01,
446            power: Some(12),
447            enabled_dampening: true,
448        };
449        settings.hard_coded_integrand = Some(IntegrandSettings::HFunctionTest(itg.clone()));
450        assert!(compare_integration(
451            &mut settings,
452            &Model::default(),
453            IntegratedPhase::Real,
454            SymComplex::new_one().into(),
455            None
456        )?);
457        Ok(())
458    }
459}
460
461#[cfg(test)]
462mod tests_inspect {
463    use crate::{
464        settings::runtime::ParameterizationMode,
465        settings::runtime::ParameterizationSettings,
466        settings::runtime::{HFunction, HFunctionSettings, SamplingSettings},
467    };
468
469    use super::*;
470
471    // Amazingly enough, a simple failing test induces a segfault on MacOS... :/
472    // #[test]
473    // fn this_test_will_not_pass() {
474    //     assert!(false);
475    // }
476
477    mod failing {
478        use super::*;
479
480        #[test]
481        fn inspect_unit_volume() -> Result<()> {
482            let mut settings = load_default_settings();
483            let mut itg = get_unit_volume_integrand();
484            itg.n_3d_momenta = 6;
485            settings.kinematics.e_cm = 1.;
486            let sampling_settings = SamplingSettings::Default(ParameterizationSettings {
487                mode: ParameterizationMode::Spherical,
488                ..Default::default()
489            });
490            settings.sampling = sampling_settings;
491
492            settings.hard_coded_integrand = Some(IntegrandSettings::UnitVolume(itg.clone()));
493            assert!(compare_inspect(
494                &mut settings,
495                &Model::default(),
496                vec![0.1, 0.2, 0.3, 0.4, 0.5, 0.6],
497                &[0],
498                true,
499                Complex::new(1.3793965770302298e3, 0.0)
500            )?);
501
502            itg.n_3d_momenta = 9;
503            settings.kinematics.e_cm = 100.;
504
505            let sampling_settings = SamplingSettings::Default(ParameterizationSettings {
506                mode: ParameterizationMode::HyperSpherical,
507                ..Default::default()
508            });
509
510            settings.sampling = sampling_settings;
511            settings.hard_coded_integrand = Some(IntegrandSettings::UnitVolume(itg.clone()));
512            assert!(compare_inspect(
513                &mut settings,
514                &Model::default(),
515                vec![0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9],
516                &[0],
517                true,
518                Complex::new(4.792927924134406e-45, 0.0)
519            )?);
520            Ok(())
521        }
522
523        #[test]
524        fn inspect_h_function_test() -> Result<()> {
525            let mut settings = load_default_settings();
526            let mut itg = get_h_function_test_integrand();
527            settings.kinematics.e_cm = 1.;
528
529            itg.h_function = HFunctionSettings {
530                function: HFunction::PolyLeftRightExponential,
531                sigma: 0.01,
532                power: Some(12),
533                enabled_dampening: true,
534            };
535            settings.hard_coded_integrand = Some(IntegrandSettings::HFunctionTest(itg.clone()));
536            assert!(compare_inspect(
537                &mut settings,
538                &Model::default(),
539                vec![0.2188450233532342,],
540                &[0],
541                false,
542                Complex::new(1.4016882047579115e-34, 0.0)
543            )?);
544
545            itg.h_function = HFunctionSettings {
546                function: HFunction::PolyLeftRightExponential,
547                sigma: 0.3,
548                power: Some(9),
549                enabled_dampening: false,
550            };
551            settings.hard_coded_integrand = Some(IntegrandSettings::HFunctionTest(itg.clone()));
552            assert!(compare_inspect(
553                &mut settings,
554                &Model::default(),
555                vec![0.2188450233532342,],
556                &[0],
557                false,
558                Complex::new(3.112977432926161e-4, 0.0)
559            )?);
560            Ok(())
561        }
562    }
563}
gammalooprs/tests.rs

gammalooprs/
tests.rs
