From 3835b1d04444b585fa97a1d201f2e649fcde7016 Mon Sep 17 00:00:00 2001 From: Mathijs Fraaije Date: Wed, 27 Nov 2024 11:28:40 +0100 Subject: [PATCH 1/6] float refactor --- Cargo.toml | 2 +- benches/gamma_bench.rs | 6 +- benches/graph_bench.rs | 6 +- src/float.rs | 407 ++++++++++++++++++++++++++++++++++++++--- src/lib.rs | 145 +++------------ src/matrix.rs | 316 ++++++++++++++++---------------- src/mimic_rng.rs | 18 +- src/preprocessing.rs | 44 ++--- src/sampling.rs | 194 +++++++++++--------- src/vector.rs | 152 ++++++++------- tests/test_prec_eq.rs | 48 +++-- tests/triangle.rs | 40 ++-- 12 files changed, 840 insertions(+), 538 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 6b29c10..8ee7fad 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -23,9 +23,9 @@ lto = "fat" smallvec = "1.13.2" ahash = "0.8.11" num = { version = "0.4.1", features = ["serde"] } -f128 = { git = "https://github.com/benruijl/f128" } statrs = "0.16.0" itertools = "0.12.1" rand = "0.8.5" criterion = "0.5.1" serde = { version = "1.0", features = ["derive"] } +ref-ops = "0.2.5" diff --git a/benches/gamma_bench.rs b/benches/gamma_bench.rs index 7a2c1f9..b93f0f8 100644 --- a/benches/gamma_bench.rs +++ b/benches/gamma_bench.rs @@ -1,13 +1,13 @@ use criterion::{criterion_group, criterion_main, Criterion}; -use momtrop::gamma::inverse_gamma_lr; +use momtrop::gamma::inverse_gamma_lr_impl; fn criterion_benchmark(c: &mut Criterion) { let mut group = c.benchmark_group("gamma sampling benchmarks"); for omega in [0.5, 1.0, 2.0, 10.0] { for p in [0.1, 0.3, 0.5, 0.7, 0.9] { - group.bench_function(&format!("omega = {omega}, p = {p}"), |b| { - b.iter(|| inverse_gamma_lr(omega, p, 50, 5.0)) + group.bench_function(format!("omega = {omega}, p = {p}"), |b| { + b.iter(|| inverse_gamma_lr_impl(omega, p, 50, 5.0)) }); } } diff --git a/benches/graph_bench.rs b/benches/graph_bench.rs index b2410f7..80baab6 100644 --- a/benches/graph_bench.rs +++ b/benches/graph_bench.rs @@ -38,7 +38,11 @@ fn criterion_benchmark(c: &mut Criterion) { let mut rng = rand::rngs::StdRng::seed_from_u64(69); let p1 = Vector::from_array([3.0, 4.0, 5.0]); let p2 = Vector::from_array([6.0, 7.0, 8.0]); - let edge_data = vec![(None, Vector::new()), (None, p1), (None, (&p1 + &p2))]; + let edge_data = vec![ + (None, Vector::new_from_num(&0.0)), + (None, p1), + (None, (&p1 + &p2)), + ]; let x_space_point = vec![rng.r#gen(); sampler.get_dimension()]; diff --git a/src/float.rs b/src/float.rs index ad2725c..d9fb8a5 100644 --- a/src/float.rs +++ b/src/float.rs @@ -1,33 +1,388 @@ -use std::fmt::{Debug, Display, LowerExp}; -use std::iter::Sum; -use std::ops::{Add, AddAssign, Div, MulAssign, Sub, SubAssign}; +use core::f64; +use ref_ops::{RefAdd, RefDiv, RefMul, RefNeg, RefSub}; +use std::fmt::Debug; +use std::ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub, SubAssign}; -use num::traits::{Float, FloatConst}; -use num::traits::{Inv, One, Zero}; - -pub trait FloatLike: - From - + AddAssign - + Div - + Zero - + One - + Inv - + PartialOrd - + PartialEq - + Display +#[allow(clippy::wrong_self_convention)] +pub trait MomTropFloat: + for<'a> RefAdd<&'a Self, Output = Self> + + for<'a> RefSub<&'a Self, Output = Self> + + for<'a> RefMul<&'a Self, Output = Self> + + for<'a> RefDiv<&'a Self, Output = Self> + + for<'a> RefAdd + + for<'a> RefSub + + for<'a> RefMul + + for<'a> RefDiv + + Add + Sub + + Mul + + Div + + for<'a> Add<&'a Self, Output = Self> + + for<'a> Sub<&'a Self, Output = Self> + + for<'a> Mul<&'a Self, Output = Self> + + for<'a> Div<&'a Self, Output = Self> + + for<'a> RefNeg + + for<'a> AddAssign<&'a Self> + + for<'a> SubAssign<&'a Self> + + for<'a> MulAssign<&'a Self> + + Neg + Clone - + SubAssign - + Float - + MulAssign - + Add - + FloatConst - + Into - + Sum + + PartialEq + Debug - + LowerExp + + PartialOrd { + fn one(&self) -> Self; + fn ln(&self) -> Self; + fn exp(&self) -> Self; + fn cos(&self) -> Self; + fn sin(&self) -> Self; + fn powf(&self, power: &Self) -> Self; + fn sqrt(&self) -> Self; + fn from_isize(&self, value: isize) -> Self; + fn from_f64(&self, value: f64) -> Self; + fn inv(&self) -> Self; + fn to_f64(&self) -> f64; + fn zero(&self) -> Self; + fn abs(&self) -> Self; + #[allow(non_snake_case)] + fn PI(&self) -> Self; +} + +impl MomTropFloat for f64 { + fn ln(&self) -> Self { + f64::ln(*self) + } + + fn exp(&self) -> Self { + f64::exp(*self) + } + + fn cos(&self) -> Self { + f64::cos(*self) + } + + fn sin(&self) -> Self { + f64::sin(*self) + } + + fn powf(&self, power: &Self) -> Self { + f64::powf(*self, *power) + } + + fn from_f64(&self, value: f64) -> Self { + value + } + + fn from_isize(&self, value: isize) -> Self { + value as f64 + } + + fn sqrt(&self) -> Self { + f64::sqrt(*self) + } + + fn inv(&self) -> Self { + 1.0 / self + } + + fn to_f64(&self) -> f64 { + *self + } + + fn PI(&self) -> Self { + f64::consts::PI + } + + fn zero(&self) -> Self { + 0.0 + } + + fn one(&self) -> Self { + 1.0 + } + + fn abs(&self) -> Self { + f64::abs(*self) + } } -impl FloatLike for f64 {} -impl FloatLike for f128::f128 {} +//#[derive(Clone, Copy)] +//pub struct F(pub T); +// +//impl Add<&F> for &F { +// type Output = F; +// +// fn add(self, rhs: &F) -> Self::Output { +// F(self.0.ref_add(&rhs.0)) +// } +//} +// +//impl Sub<&F> for &F { +// type Output = F; +// +// fn sub(self, rhs: &F) -> Self::Output { +// F(self.0.ref_sub(&rhs.0)) +// } +//} +// +//impl Mul<&F> for &F { +// type Output = F; +// +// fn mul(self, rhs: &F) -> Self::Output { +// F(self.0.ref_mul(&rhs.0)) +// } +//} +// +//impl Div<&F> for &F { +// type Output = F; +// +// fn div(self, rhs: &F) -> Self::Output { +// println!("calling ref_div"); +// F(self.0.ref_div(&rhs.0)) +// } +//} +// +//// additional addition traits +//impl Add> for &F { +// type Output = F; +// +// fn add(self, rhs: F) -> Self::Output { +// self + &rhs +// } +//} +// +//impl Add<&F> for F { +// type Output = F; +// +// fn add(self, rhs: &F) -> Self::Output { +// &self + rhs +// } +//} +// +//impl Add> for F { +// type Output = F; +// +// fn add(self, rhs: F) -> Self::Output { +// &self + &rhs +// } +//} +// +//// additional subtraction traits +//impl Sub> for &F { +// type Output = F; +// +// fn sub(self, rhs: F) -> Self::Output { +// println!("recursing"); +// self - &rhs +// } +//} +// +//impl Sub<&F> for F { +// type Output = F; +// +// fn sub(self, rhs: &F) -> Self::Output { +// println!("recursing"); +// &self - rhs +// } +//} +// +//impl Sub> for F { +// type Output = F; +// +// fn sub(self, rhs: F) -> Self::Output { +// println!("recursing"); +// &self - &rhs +// } +//} +// +//// additional multiplication traits +//impl Mul> for &F { +// type Output = F; +// +// fn mul(self, rhs: F) -> Self::Output { +// println!("recursing"); +// self * &rhs +// } +//} +// +//impl Mul<&F> for F { +// type Output = F; +// +// fn mul(self, rhs: &F) -> Self::Output { +// println!("recursing"); +// &self * rhs +// } +//} +// +//impl Mul> for F { +// type Output = F; +// +// fn mul(self, rhs: F) -> Self::Output { +// println!("recursing"); +// &self * &rhs +// } +//} +// +//impl Div> for &F { +// type Output = F; +// +// fn div(self, rhs: F) -> Self::Output { +// println!("no_ref ref"); +// self / &rhs +// } +//} +// +//impl Div<&F> for F { +// type Output = F; +// +// fn div(self, rhs: &F) -> Self::Output { +// println!("ref no_ref"); +// &self / rhs +// } +//} +// +//impl Div> for F { +// type Output = F; +// +// fn div(self, rhs: F) -> Self::Output { +// println!("no_ref no_ref"); +// &self / &rhs +// } +//} +// +//impl Neg for &F { +// type Output = F; +// +// fn neg(self) -> Self::Output { +// self.ref_neg() +// } +//} +// +//impl Neg for F { +// type Output = F; +// +// fn neg(self) -> Self::Output { +// -&self +// } +//} +// +//impl AddAssign<&F> for F { +// fn add_assign(&mut self, rhs: &F) { +// *self = &*self + rhs; +// } +//} +// +//impl AddAssign> for F { +// fn add_assign(&mut self, rhs: F) { +// *self = &*self + rhs; +// } +//} +// +//impl SubAssign<&F> for F { +// fn sub_assign(&mut self, rhs: &F) { +// println!("recursing"); +// *self = &*self - rhs; +// } +//} +// +//impl SubAssign> for F { +// fn sub_assign(&mut self, rhs: F) { +// println!("recursing"); +// *self = &*self - rhs; +// } +//} +// +//impl MulAssign<&F> for F { +// fn mul_assign(&mut self, rhs: &F) { +// *self = &*self * rhs +// } +//} +// +//impl MulAssign> for F { +// fn mul_assign(&mut self, rhs: F) { +// *self = &*self * rhs +// } +//} +// +//impl PartialEq for F { +// fn eq(&self, other: &Self) -> bool { +// self.0 == other.0 +// } +//} +// +//impl PartialOrd for F { +// fn partial_cmp(&self, other: &Self) -> Option { +// self.0.partial_cmp(&other.0) +// } +//} +// +//impl Debug for F { +// fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { +// self.0.fmt(f) +// } +//} +// +//impl F { +// pub fn zero(&self) -> Self { +// Self(self.0.from_isize(0)) +// } +// +// pub fn one(&self) -> Self { +// Self(self.0.from_isize(1)) +// } +// +// pub fn from_isize(&self, value: isize) -> Self { +// Self(self.0.from_isize(value)) +// } +// +// pub fn from_f64(&self, value: f64) -> Self { +// Self(self.0.from_f64(value)) +// } +// +// pub fn abs(&self) -> Self { +// if self < &self.zero() { +// -self +// } else { +// self.clone() +// } +// } +// +// pub fn sqrt(&self) -> Self { +// Self(self.0.sqrt()) +// } +// +// pub fn ln(&self) -> Self { +// Self(self.0.ln()) +// } +// +// pub fn inv(&self) -> Self { +// Self(self.0.inv()) +// } +// +// pub fn to_f64(&self) -> f64 { +// self.0.to_f64() +// } +// +// #[allow(non_snake_case)] +// pub fn PI(&self) -> Self { +// Self(self.0.PI()) +// } +// +// pub fn cos(&self) -> Self { +// Self(self.0.cos()) +// } +// +// pub fn sin(&self) -> Self { +// Self(self.0.sin()) +// } +// +// pub fn powf(&self, power: &Self) -> Self { +// Self(self.0.powf(&power.0)) +// } +//} +// +//impl From for F { +// fn from(value: T) -> Self { +// Self(value) +// } +//} diff --git a/src/lib.rs b/src/lib.rs index 9813c52..d56656d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -3,8 +3,7 @@ use std::iter::repeat_with; #[cfg(feature = "log")] use log::Logger; -use f128::f128; -use float::FloatLike; +use float::MomTropFloat; use itertools::Itertools; use matrix::{DecompositionResult, SquareMatrix}; use preprocessing::{TropicalGraph, TropicalSubgraphTable}; @@ -13,7 +12,7 @@ use sampling::{sample, SamplingError}; use serde::{Deserialize, Serialize}; use vector::Vector; -mod float; +pub mod float; pub mod gamma; #[cfg(feature = "log")] pub mod log; @@ -45,13 +44,12 @@ impl Default for TropicalSamplingSettings { } } -#[cfg(test)] -fn assert_approx_eq(res: f64, target: f64, tolerance: f64) { +pub fn assert_approx_eq(res: &T, target: &T, tolerance: &T) { if approx_eq(res, target, tolerance) { } else { panic!( - "assert_approx_eq failed: \n{:+e} != \n{:+e} with tolerance {:+e}", - res, target, tolerance + "assert_approx_eq failed: \n{:?} != \n{:?} with tolerance {:?}", + &res, &target, &tolerance ) } } @@ -68,7 +66,7 @@ pub struct Edge { } #[derive(Debug, Clone)] -pub struct TropicalSampleResult { +pub struct TropicalSampleResult { pub loop_momenta: Vec>, pub u_trop: T, pub v_trop: T, @@ -79,7 +77,7 @@ pub struct TropicalSampleResult { } #[derive(Debug, Clone)] -pub struct Metadata { +pub struct Metadata { pub q_vectors: Vec>, pub lambda: T, pub l_matrix: SquareMatrix, @@ -88,33 +86,6 @@ pub struct Metadata { pub shift: Vec>, } -impl Metadata { - pub fn downcast(&self) -> Metadata { - Metadata { - q_vectors: self.q_vectors.iter().map(|v| v.downcast()).collect_vec(), - lambda: self.lambda.into(), - l_matrix: self.l_matrix.downcast(), - decompoisiton_result: self.decompoisiton_result.downcast(), - u_vectors: self.u_vectors.iter().map(|v| v.downcast()).collect_vec(), - shift: self.shift.iter().map(|v| v.downcast()).collect(), - } - } -} - -impl TropicalSampleResult { - pub fn downcast(&self) -> TropicalSampleResult { - TropicalSampleResult { - loop_momenta: self.loop_momenta.iter().map(|v| v.downcast()).collect_vec(), - u_trop: self.u_trop.into(), - v_trop: self.v_trop.into(), - u: self.u.into(), - v: self.v.into(), - jacobian: self.jacobian.into(), - metadata: self.metadata.as_ref().map(|metadata| metadata.downcast()), - } - } -} - impl Graph { pub fn build_sampler( self, @@ -131,12 +102,11 @@ impl Graph { } } -#[cfg(test)] -fn approx_eq(res: f64, target: f64, tolerance: f64) -> bool { - if target == 0.0 { - res.abs() < tolerance +fn approx_eq(res: &T, target: &T, tolerance: &T) -> bool { + if target == &res.zero() { + &res.abs() < tolerance } else { - ((res - target) / target).abs() < tolerance + &((res.ref_sub(target)) / target).abs() < tolerance } } @@ -147,14 +117,17 @@ pub struct SampleGenerator { } impl SampleGenerator { - pub fn generate_sample_from_x_space_point<#[cfg(feature = "log")] L: Logger>( + pub fn generate_sample_from_x_space_point< + T: MomTropFloat, + #[cfg(feature = "log")] L: Logger, + >( &self, - x_space_point: &[f64], - edge_data: Vec<(Option, vector::Vector)>, + x_space_point: &[T], + edge_data: Vec<(Option, vector::Vector)>, settings: &TropicalSamplingSettings, #[cfg(feature = "log")] logger: &L, - ) -> Result, SamplingError> { - let sample = sample( + ) -> Result, SamplingError> { + sample( &self.table, x_space_point, &self.loop_signature, @@ -162,90 +135,24 @@ impl SampleGenerator { settings, #[cfg(feature = "log")] logger, - ); - - match sample { - Ok(sample) => Ok(sample), - Err(sampling_error) => { - if settings.upcast_on_failure { - match sampling_error { - SamplingError::MatrixError(_matrix_error) => { - let res = self.generate_sample_f128_from_x_space_point( - x_space_point, - edge_data, - settings, - #[cfg(feature = "log")] - logger, - ); - - res.map(|res| res.downcast()) - } - } - } else { - Err(sampling_error) - } - } - } - } - - pub fn generate_sample_from_rng( - &self, - edge_data: Vec<(Option, vector::Vector)>, - settings: &TropicalSamplingSettings, - rng: &mut R, - #[cfg(feature = "log")] logger: &L, - ) -> Result, SamplingError> { - let num_vars = self.get_dimension(); - let x_space_point = repeat_with(|| rng.gen::()) - .take(num_vars) - .collect_vec(); - - self.generate_sample_from_x_space_point( - &x_space_point, - edge_data, - settings, - #[cfg(feature = "log")] - logger, - ) - } - - pub fn generate_sample_f128_from_x_space_point<#[cfg(feature = "log")] L: Logger>( - &self, - x_space_point: &[f64], - edge_data: Vec<(Option, vector::Vector)>, - settings: &TropicalSamplingSettings, - #[cfg(feature = "log")] logger: &L, - ) -> Result, SamplingError> { - let upcasted_xspace_point = x_space_point.iter().copied().map(f128::new).collect_vec(); - let upcasted_edge_data = edge_data - .into_iter() - .map(|(mass, edge_shift)| (mass.map(f128::new), edge_shift.upcast())) - .collect_vec(); - - sample( - &self.table, - &upcasted_xspace_point, - &self.loop_signature, - &upcasted_edge_data, - settings, - #[cfg(feature = "log")] - logger, ) } - pub fn generate_sample_f128_from_rng( + pub fn generate_sample_from_rng( &self, - edge_data: Vec<(Option, vector::Vector)>, + edge_data: Vec<(Option, vector::Vector)>, settings: &TropicalSamplingSettings, rng: &mut R, #[cfg(feature = "log")] logger: &L, - ) -> Result, SamplingError> { + ) -> Result, SamplingError> { + let const_builder = edge_data[0].1.zero(); + let num_vars = self.get_dimension(); - let x_space_point = repeat_with(|| rng.gen::()) + let x_space_point = repeat_with(|| const_builder.from_f64(rng.gen::())) .take(num_vars) .collect_vec(); - self.generate_sample_f128_from_x_space_point( + self.generate_sample_from_x_space_point( &x_space_point, edge_data, settings, diff --git a/src/matrix.rs b/src/matrix.rs index 17ca17a..9a49bd4 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -1,8 +1,7 @@ use smallvec::SmallVec; use std::ops::{Add, Index, IndexMut, Mul, Sub}; -use crate::{float::FloatLike, TropicalSamplingSettings}; -use f128::f128; +use crate::{float::MomTropFloat, TropicalSamplingSettings}; #[derive(Debug, Clone, Default)] /// square symmetric matrix for use in the tropical sampling algorithm @@ -31,16 +30,18 @@ impl IndexMut<(usize, usize)> for SquareMatrix { } } -impl<'a, T: FloatLike> Mul<&'a SquareMatrix> for &'a SquareMatrix { +// silly clippy +#[allow(clippy::suspicious_arithmetic_impl)] +impl<'a, T: MomTropFloat> Mul<&'a SquareMatrix> for &'a SquareMatrix { type Output = SquareMatrix; fn mul(self, rhs: Self) -> Self::Output { - let mut result = SquareMatrix::new_zeros(self.dim); + let mut result = self.new_zeros(self.dim); for row in 0..self.dim { for col in 0..self.dim { for k in 0..self.dim { - result[(row, col)] += self[(row, k)] * rhs[(k, col)]; + result[(row, col)] += &self[(row, k)].ref_mul(&rhs[(k, col)]); } } } @@ -49,15 +50,15 @@ impl<'a, T: FloatLike> Mul<&'a SquareMatrix> for &'a SquareMatrix { } } -impl Mul for SquareMatrix { +impl Mul for SquareMatrix { type Output = SquareMatrix; fn mul(self, rhs: T) -> Self::Output { - let mut result = SquareMatrix::new_zeros(self.dim); + let mut result = self.new_zeros(self.dim); for row in 0..self.dim { for col in 0..self.dim { - result[(row, col)] = self[(row, col)] * rhs; + result[(row, col)] = self[(row, col)].ref_mul(&rhs); } } @@ -65,15 +66,15 @@ impl Mul for SquareMatrix { } } -impl<'a, T: FloatLike> Add<&'a SquareMatrix> for &'a SquareMatrix { +impl<'a, T: MomTropFloat> Add<&'a SquareMatrix> for &'a SquareMatrix { type Output = SquareMatrix; fn add(self, rhs: Self) -> Self::Output { - let mut result = SquareMatrix::new_zeros(self.dim); + let mut result = self.new_zeros(self.dim); for row in 0..self.dim { for col in 0..self.dim { - result[(row, col)] = self[(row, col)] + rhs[(row, col)]; + result[(row, col)] = self[(row, col)].ref_add(&rhs[(row, col)]); } } @@ -81,15 +82,15 @@ impl<'a, T: FloatLike> Add<&'a SquareMatrix> for &'a SquareMatrix { } } -impl<'a, T: FloatLike> Sub<&'a SquareMatrix> for &'a SquareMatrix { +impl<'a, T: MomTropFloat> Sub<&'a SquareMatrix> for &'a SquareMatrix { type Output = SquareMatrix; fn sub(self, rhs: Self) -> Self::Output { - let mut result = SquareMatrix::new_zeros(self.dim); + let mut result = self.new_zeros(self.dim); for row in 0..self.dim { for col in 0..self.dim { - result[(row, col)] = self[(row, col)] - rhs[(row, col)]; + result[(row, col)] = self[(row, col)].ref_sub(&rhs[(row, col)]); } } @@ -97,14 +98,22 @@ impl<'a, T: FloatLike> Sub<&'a SquareMatrix> for &'a SquareMatrix { } } -impl SquareMatrix { +impl SquareMatrix { #[must_use] - pub fn new_zeros(dim: usize) -> Self { + pub fn new_zeros(&self, dim: usize) -> Self { Self { - data: SmallVec::from_elem(T::zero(), dim * dim), + data: SmallVec::from_elem(self.data[0].zero(), dim * dim), dim, } } + + pub fn new_zeros_from_num(builder: &T, dim: usize) -> Self { + Self { + data: SmallVec::from_elem(builder.zero(), dim * dim), + dim, + } + } + /// Performs operations on a matrix for tropical sampling /// # Errors /// Returns an error if the matrix is not invertible @@ -113,54 +122,55 @@ impl SquareMatrix { settings: &TropicalSamplingSettings, ) -> Result, MatrixError> { if settings.print_debug_info { - println!("starting cholesky"); + println!("starting cholesky") } + let const_builder = &self.data[0]; // start cholesky decomposition - let mut q: SquareMatrix = SquareMatrix::new_zeros(self.dim); + let mut q = self.new_zeros(self.dim); for i in 0..self.dim { - let mut diagonal_entry_squared = self[(i, i)]; + let mut diagonal_entry_squared = self[(i, i)].clone(); for j in 0..i { - diagonal_entry_squared -= q[(i, j)] * q[(i, j)]; + diagonal_entry_squared -= &q[(i, j)].ref_mul(&q[(i, j)]); } let diagonal_entry = diagonal_entry_squared.sqrt(); - q[(i, i)] = diagonal_entry; + q[(i, i)] = diagonal_entry.clone(); for j in i + 1..self.dim { - let mut entry = self[(i, j)]; + let mut entry = self[(i, j)].clone(); for k in 0..i { - entry -= q[(i, k)] * q[(j, k)]; + entry -= &q[(i, k)].ref_mul(&q[(j, k)]); } - q[(j, i)] = entry / diagonal_entry; + q[(j, i)] = entry / &diagonal_entry; } } // end cholesky decomposition // compute the determinant of Q and store the invesrses of the diagonal elements - let mut det_q = T::one(); + let mut det_q = const_builder.one(); let mut inverse_diagonal_entries = SmallVec::<[T; 6]>::new(); for i in 0..self.dim { - let q_ii = q[(i, i)]; - det_q *= q_ii; + let q_ii = q[(i, i)].clone(); + det_q *= &q_ii; inverse_diagonal_entries.push(q_ii.inv()); } - (0..self.dim).fold(T::one(), |acc, i| acc * q[(i, i)]); - let determinant = det_q * det_q; + (0..self.dim).fold(const_builder.one(), |acc, i| acc * &q[(i, i)]); + let determinant = det_q.ref_mul(&det_q); - if det_q == T::zero() { + if det_q == const_builder.zero() { return Err(MatrixError::ZeroDet); } // the matrix N is defined through Q = D(I + N) - let mut n_matrix: SquareMatrix = SquareMatrix::new_zeros(self.dim); + let mut n_matrix = self.new_zeros(self.dim); for row in 1..self.dim { - let inverse_diagonal_element = inverse_diagonal_entries[row]; + let inverse_diagonal_element = &inverse_diagonal_entries[row]; for col in 0..row { - n_matrix[(row, col)] = inverse_diagonal_element * q[(row, col)]; + n_matrix[(row, col)] = inverse_diagonal_element.ref_mul(&q[(row, col)]); } } @@ -180,31 +190,32 @@ impl SquareMatrix { powers_of_n.push(last_power_of_n * first_power_of_n); } - let n_sum = powers_of_n.iter().enumerate().fold( - SquareMatrix::new_zeros(self.dim), - |acc, (i, mat)| { - if i % 2 == 0 { - &acc - mat - } else { - &acc + mat - } - }, - ); + let n_sum = + powers_of_n + .iter() + .enumerate() + .fold(self.new_zeros(self.dim), |acc, (i, mat)| { + if i % 2 == 0 { + &acc - mat + } else { + &acc + mat + } + }); let mut inverse_q = n_sum; for row in 0..self.dim { - inverse_q[(row, row)] += T::one(); + inverse_q[(row, row)] += &const_builder.one(); for col in 0..self.dim { - inverse_q[(row, col)] *= inverse_diagonal_entries[col]; + inverse_q[(row, col)] *= &inverse_diagonal_entries[col]; } } - let mut q_transposed_inverse = SquareMatrix::new_zeros(self.dim); - let mut q_transposed = SquareMatrix::new_zeros(self.dim); + let mut q_transposed_inverse = self.new_zeros(self.dim); + let mut q_transposed = self.new_zeros(self.dim); for row in 0..self.dim { for col in 0..self.dim { - q_transposed_inverse[(row, col)] = inverse_q[(col, row)]; - q_transposed[(row, col)] = q[(col, row)]; + q_transposed_inverse[(row, col)] = inverse_q[(col, row)].clone(); + q_transposed[(row, col)] = q[(col, row)].clone(); } } @@ -216,15 +227,15 @@ impl SquareMatrix { } let approx_idendity = &inverse * self; - let true_identity = Self::new_identity(self.dim); + let true_identity = self.new_identity(self.dim); let zero = &approx_idendity - &true_identity; let error = zero.l21_norm(); if settings.print_debug_info { - println!("error: {:+16e}", error); + println!("error: {:?}", error); } - if Into::::into(error) > tolerance { + if error > error.from_f64(tolerance) { if settings.print_debug_info { println!("Inversion unstable"); } @@ -241,27 +252,31 @@ impl SquareMatrix { }) } - fn new_identity(dim: usize) -> Self { - let mut res = Self::new_zeros(dim); + fn new_identity(&self, dim: usize) -> Self { + let mut res = self.new_zeros(dim); for i in 0..dim { - res[(i, i)] = T::one(); + res[(i, i)] = self.data[0].one(); } res } fn l21_norm(&self) -> T { - let mut res = T::zero(); + let mut res = self.data[0].zero(); for j in 0..self.dim { - let mut vec_norm = T::zero(); + let mut vec_norm = res.zero(); for i in 0..self.dim { - vec_norm += self[(i, j)] * self[(i, j)]; + vec_norm += &self[(i, j)].ref_mul(&self[(i, j)]); } - res += vec_norm.sqrt(); + res += &vec_norm.sqrt(); } res } + + pub fn zero(&self) -> T { + self.data[0].zero() + } } impl SquareMatrix { @@ -274,34 +289,14 @@ impl SquareMatrix { } } -impl SquareMatrix { - pub fn downcast(&self) -> SquareMatrix { - SquareMatrix { - data: self.data.iter().map(|&x| x.into()).collect(), - dim: self.dim, - } - } -} - #[derive(Debug, Clone)] -pub struct DecompositionResult { +pub struct DecompositionResult { pub determinant: T, pub inverse: SquareMatrix, pub q_transposed: SquareMatrix, pub q_transposed_inverse: SquareMatrix, } -impl DecompositionResult { - pub fn downcast(&self) -> DecompositionResult { - DecompositionResult { - determinant: self.determinant.into(), - q_transposed: self.q_transposed.downcast(), - inverse: self.inverse.downcast(), - q_transposed_inverse: self.q_transposed_inverse.downcast(), - } - } -} - #[cfg(test)] mod tests { use super::*; @@ -309,10 +304,23 @@ mod tests { const EPSILON: f64 = 1e-12; + fn builder_matrix_f64() -> SquareMatrix { + SquareMatrix { + data: smallvec::SmallVec::from_elem(1.0, 1), + dim: 1, + } + } + #[test] fn test_decompose_for_tropical_2x2() { - let mut test_matrix = SquareMatrix::new_zeros(2); - let settings = TropicalSamplingSettings::default(); + println!("starting test"); + let mut test_matrix = builder_matrix_f64().new_zeros(2); + println!("build matrix"); + + let settings = TropicalSamplingSettings { + print_debug_info: true, + ..Default::default() + }; test_matrix[(0, 0)] = 2.0; test_matrix[(1, 1)] = 4.0; @@ -321,45 +329,45 @@ mod tests { let decomposition_result = test_matrix.decompose_for_tropical(&settings).unwrap(); - assert_approx_eq(decomposition_result.determinant, 7.0, EPSILON); + assert_approx_eq(&decomposition_result.determinant, &7.0, &EPSILON); let identity = &test_matrix * &decomposition_result.inverse; - assert_approx_eq(identity[(0, 0)], 1.0, EPSILON); - assert_approx_eq(identity[(0, 1)], 0.0, EPSILON); - assert_approx_eq(identity[(1, 0)], 0.0, EPSILON); - assert_approx_eq(identity[(1, 1)], 1.0, EPSILON); + assert_approx_eq(&identity[(0, 0)], &1.0, &EPSILON); + assert_approx_eq(&identity[(0, 1)], &0.0, &EPSILON); + assert_approx_eq(&identity[(1, 0)], &0.0, &EPSILON); + assert_approx_eq(&identity[(1, 1)], &1.0, &EPSILON); } - #[test] - fn test_decompose_for_tropical_f128_2x2() { - let mut test_matrix = SquareMatrix::new_zeros(2); - let settings = TropicalSamplingSettings::default(); + // #[test] + // fn test_decompose_for_tropical_f128_2x2() { + // let mut test_matrix = SquareMatrix::new_zeros(2); + // let settings = TropicalSamplingSettings::default(); - test_matrix[(0, 0)] = f128::new(2.0); - test_matrix[(1, 1)] = f128::new(4.0); - test_matrix[(0, 1)] = f128::new(1.0); - test_matrix[(1, 0)] = f128::new(1.0); + // test_matrix[(0, 0)] = f128::new(2.0); + // test_matrix[(1, 1)] = f128::new(4.0); + // test_matrix[(0, 1)] = f128::new(1.0); + // test_matrix[(1, 0)] = f128::new(1.0); - let decomposition_result = test_matrix.decompose_for_tropical(&settings).unwrap(); + // let decomposition_result = test_matrix.decompose_for_tropical(&settings).unwrap(); - assert_approx_eq( - Into::::into(decomposition_result.determinant), - 7.0, - EPSILON, - ); + // assert_approx_eq( + // Into::::into(decomposition_result.determinant), + // 7.0, + // EPSILON, + // ); - let identity = &test_matrix * &decomposition_result.inverse; + // let identity = &test_matrix * &decomposition_result.inverse; - assert_approx_eq(Into::::into(identity[(0, 0)]), 1.0, EPSILON); - assert_approx_eq(Into::::into(identity[(0, 1)]), 0.0, EPSILON); - assert_approx_eq(Into::::into(identity[(1, 0)]), 0.0, EPSILON); - assert_approx_eq(Into::::into(identity[(1, 1)]), 1.0, EPSILON); - } + // assert_approx_eq(Into::::into(identity[(0, 0)]), 1.0, EPSILON); + // assert_approx_eq(Into::::into(identity[(0, 1)]), 0.0, EPSILON); + // assert_approx_eq(Into::::into(identity[(1, 0)]), 0.0, EPSILON); + // assert_approx_eq(Into::::into(identity[(1, 1)]), 1.0, EPSILON); + // } #[test] fn test_decompose_for_tropical_4x4() { - let mut wilson_matrix = SquareMatrix::new_zeros(4); + let mut wilson_matrix = builder_matrix_f64().new_zeros(4); let settings = TropicalSamplingSettings::default(); wilson_matrix[(0, 0)] = 5.0; @@ -386,55 +394,55 @@ mod tests { for row in 0..4 { for col in 0..4 { assert_approx_eq( - identity[(row, col)], - if row == col { 1.0 } else { 0.0 }, - EPSILON, + &identity[(row, col)], + if row == col { &1.0 } else { &0.0 }, + &EPSILON, ); } } - assert_approx_eq(decomposition_result.determinant, 1.0, EPSILON); + assert_approx_eq(&decomposition_result.determinant, &1.0, &EPSILON); } - #[test] - fn test_decompose_for_tropical_4x4_f128() { - let mut wilson_matrix = SquareMatrix::new_zeros(4); - let settings = TropicalSamplingSettings::default(); - - wilson_matrix[(0, 0)] = f128::new(5.0); - wilson_matrix[(1, 1)] = f128::new(10.); - wilson_matrix[(2, 2)] = f128::new(10.); - wilson_matrix[(3, 3)] = f128::new(10.); - wilson_matrix[(0, 1)] = f128::new(7.0); - wilson_matrix[(1, 0)] = f128::new(7.0); - wilson_matrix[(0, 2)] = f128::new(6.0); - wilson_matrix[(2, 0)] = f128::new(6.0); - wilson_matrix[(0, 3)] = f128::new(5.0); - wilson_matrix[(3, 0)] = f128::new(5.0); - wilson_matrix[(1, 2)] = f128::new(8.0); - wilson_matrix[(2, 1)] = f128::new(8.0); - wilson_matrix[(1, 3)] = f128::new(7.0); - wilson_matrix[(3, 1)] = f128::new(7.0); - wilson_matrix[(2, 3)] = f128::new(9.0); - wilson_matrix[(3, 2)] = f128::new(9.0); - - let decomposition_result = wilson_matrix.decompose_for_tropical(&settings).unwrap(); - let identity = &wilson_matrix * &decomposition_result.inverse; - - for row in 0..4 { - for col in 0..4 { - assert_approx_eq( - Into::::into(identity[(row, col)]), - if row == col { 1.0 } else { 0.0 }, - EPSILON, - ); - } - } - - assert_approx_eq( - Into::::into(decomposition_result.determinant), - 1.0, - EPSILON, - ); - } + // #[test] + // fn test_decompose_for_tropical_4x4_f128() { + // let mut wilson_matrix = SquareMatrix::new_zeros(4); + // let settings = TropicalSamplingSettings::default(); + + // wilson_matrix[(0, 0)] = f128::new(5.0); + // wilson_matrix[(1, 1)] = f128::new(10.); + // wilson_matrix[(2, 2)] = f128::new(10.); + // wilson_matrix[(3, 3)] = f128::new(10.); + // wilson_matrix[(0, 1)] = f128::new(7.0); + // wilson_matrix[(1, 0)] = f128::new(7.0); + // wilson_matrix[(0, 2)] = f128::new(6.0); + // wilson_matrix[(2, 0)] = f128::new(6.0); + // wilson_matrix[(0, 3)] = f128::new(5.0); + // wilson_matrix[(3, 0)] = f128::new(5.0); + // wilson_matrix[(1, 2)] = f128::new(8.0); + // wilson_matrix[(2, 1)] = f128::new(8.0); + // wilson_matrix[(1, 3)] = f128::new(7.0); + // wilson_matrix[(3, 1)] = f128::new(7.0); + // wilson_matrix[(2, 3)] = f128::new(9.0); + // wilson_matrix[(3, 2)] = f128::new(9.0); + + // let decomposition_result = wilson_matrix.decompose_for_tropical(&settings).unwrap(); + // let identity = &wilson_matrix * &decomposition_result.inverse; + + // for row in 0..4 { + // for col in 0..4 { + // assert_approx_eq( + // Into::::into(identity[(row, col)]), + // if row == col { 1.0 } else { 0.0 }, + // EPSILON, + // ); + // } + // } + + // assert_approx_eq( + // Into::::into(decomposition_result.determinant), + // 1.0, + // EPSILON, + // ); + // } } diff --git a/src/mimic_rng.rs b/src/mimic_rng.rs index 544d2fe..da596e7 100644 --- a/src/mimic_rng.rs +++ b/src/mimic_rng.rs @@ -1,3 +1,5 @@ +use crate::float::MomTropFloat; + /// Fake random number generator, allows you to specifiy your own random variables pub struct MimicRng<'a, T> { cache: &'a [T], @@ -5,7 +7,7 @@ pub struct MimicRng<'a, T> { tokens: Vec<&'a str>, } -impl<'a, T: Copy> MimicRng<'a, T> { +impl<'a, T> MimicRng<'a, T> { #[inline] pub fn new(cache: &'a [T]) -> Self { Self { @@ -16,8 +18,8 @@ impl<'a, T: Copy> MimicRng<'a, T> { } #[inline] - pub fn get_random_number(&mut self, debug_token: Option<&'a str>) -> T { - let random_number = self.cache[self.counter]; + pub fn get_random_number(&mut self, debug_token: Option<&'a str>) -> &'a T { + let random_number = &self.cache[self.counter]; self.counter += 1; if let Some(token) = debug_token { self.tokens.push(token); @@ -25,3 +27,13 @@ impl<'a, T: Copy> MimicRng<'a, T> { random_number } } + +impl<'a, T: MomTropFloat> MimicRng<'a, T> { + pub fn zero(&self) -> T { + self.cache[0].zero() + } + + pub fn one(&self) -> T { + self.cache[0].one() + } +} diff --git a/src/preprocessing.rs b/src/preprocessing.rs index 2d99e8d..8af5f26 100644 --- a/src/preprocessing.rs +++ b/src/preprocessing.rs @@ -6,7 +6,7 @@ use num::Zero; use serde::{Deserialize, Serialize}; use statrs::function::gamma::gamma; -use crate::{float::FloatLike, Graph, MAX_EDGES}; +use crate::{float::MomTropFloat, Graph, MAX_EDGES}; #[derive(Debug, Clone, Serialize, Deserialize)] pub struct TropicalGraph { @@ -460,27 +460,27 @@ impl TropicalSubgraphTable { /// sample an edge from a subgraph, according to the relevant probability distribution, returns the edge and the subgraph without the edge for later use /// Panics if the uniform random number is greater than one or the probability distribution is incorrectly normalized. - pub fn sample_edge( + pub fn sample_edge( &self, - uniform: T, + uniform: &T, subgraph: &TropicalSubGraphId, ) -> (usize, TropicalSubGraphId) { let edges_in_subgraph = subgraph.contains_edges(); - let j = Into::::into(self.table[subgraph.id].j_function); + let j = uniform.from_f64(self.table[subgraph.id].j_function); - let mut cum_sum = T::zero(); + let mut cum_sum = uniform.zero(); for edge in edges_in_subgraph { let graph_without_edge = subgraph.pop_edge(edge); - let p_e = Into::::into(self.table[graph_without_edge.id].j_function) - / Into::::into(j) - / Into::::into(self.table[graph_without_edge.id].generalized_dod); - cum_sum += p_e; - if cum_sum >= uniform { + let p_e = uniform.from_f64(self.table[graph_without_edge.id].j_function) + / &j + / uniform.from_f64(self.table[graph_without_edge.id].generalized_dod); + cum_sum += &p_e; + if &cum_sum >= uniform { return (edge, graph_without_edge); } } - panic!("Sampling could not sample edge, with uniform_random_number: {}, cumulative sum evaluated to: {}", uniform, cum_sum); + panic!("Sampling could not sample edge, with uniform_random_number: {:?}, cumulative sum evaluated to: {:?}", uniform, cum_sum); } pub fn get_num_variables(&self) -> usize { @@ -984,15 +984,15 @@ mod tests { let table = subgraph_table.table[i]; assert!(table.mass_momentum_spanning); assert_eq!(table.loop_number, 0); - assert_approx_eq(table.generalized_dod, 0.84, TOLERANCE); - assert_approx_eq(table.j_function, 3.030_303_030_303_03, TOLERANCE); + assert_approx_eq(&table.generalized_dod, &0.84, &TOLERANCE); + assert_approx_eq(&table.j_function, &3.030_303_030_303_03, &TOLERANCE); } let final_table = subgraph_table.table[7]; assert!(final_table.mass_momentum_spanning); assert_eq!(final_table.loop_number, 1); - assert_approx_eq(final_table.generalized_dod, 0.0, TOLERANCE); - assert_approx_eq(final_table.j_function, 10.822_510_822_510_82, TOLERANCE); + assert_approx_eq(&final_table.generalized_dod, &0.0, &TOLERANCE); + assert_approx_eq(&final_table.j_function, &10.822_510_822_510_82, &TOLERANCE); } #[test] @@ -1048,23 +1048,23 @@ mod tests { let table = subgraph_table.table[i]; assert!(table.mass_momentum_spanning); assert_eq!(table.loop_number, 0); - assert_approx_eq(table.generalized_dod, 1.68, TOLERANCE); - assert_approx_eq(table.j_function, 1.0, TOLERANCE); + assert_approx_eq(&table.generalized_dod, &1.68, &TOLERANCE); + assert_approx_eq(&table.j_function, &1.0, &TOLERANCE); } for i in [3, 5, 6] { let table = subgraph_table.table[i]; assert!(table.mass_momentum_spanning); assert_eq!(table.loop_number, 1); - assert_approx_eq(table.generalized_dod, 0.84, TOLERANCE); - assert_approx_eq(table.j_function, 1.190_476_190_476_19, TOLERANCE); + assert_approx_eq(&table.generalized_dod, &0.84, &TOLERANCE); + assert_approx_eq(&table.j_function, &1.190_476_190_476_19, &TOLERANCE); } let final_table = subgraph_table.table[7]; assert!(final_table.mass_momentum_spanning); assert_eq!(final_table.loop_number, 2); - assert_approx_eq(final_table.generalized_dod, 0.0, TOLERANCE); - assert_approx_eq(final_table.j_function, 4.251_700_680_272_108, TOLERANCE); + assert_approx_eq(&final_table.generalized_dod, &0.0, &TOLERANCE); + assert_approx_eq(&final_table.j_function, &4.251_700_680_272_108, &TOLERANCE); } #[test] @@ -1137,6 +1137,6 @@ mod tests { let i_tr = subgraph_table.table.last().unwrap().j_function; - assert_approx_eq(i_tr, 1_818.303_855_640_347_1, TOLERANCE); + assert_approx_eq(&i_tr, &1_818.303_855_640_347_1, &TOLERANCE); } } diff --git a/src/sampling.rs b/src/sampling.rs index 71c3ec7..fd6fb76 100644 --- a/src/sampling.rs +++ b/src/sampling.rs @@ -6,13 +6,13 @@ use crate::log::Logger; use crate::matrix::{MatrixError, SquareMatrix}; use crate::mimic_rng::MimicRng; use crate::vector::Vector; -use crate::{float::FloatLike, gamma::inverse_gamma_lr}; +use crate::{float::MomTropFloat, gamma::inverse_gamma_lr}; use crate::{Metadata, TropicalSampleResult, TropicalSamplingSettings}; -fn box_muller(x1: T, x2: T) -> (T, T) { - let r = (-Into::::into(2.) * x1.ln()).sqrt(); - let theta = Into::::into(2.) * T::PI() * x2; - (r * theta.cos(), r * theta.sin()) +fn box_muller(x1: &T, x2: &T) -> (T, T) { + let r = (-x1.from_isize(2) * x1.ln()).sqrt(); + let theta = x1.from_isize(2) * x1.PI() * x2; + (theta.cos() * &r, theta.sin() * &r) } #[derive(Debug, Clone, Copy)] @@ -20,7 +20,7 @@ pub enum SamplingError { MatrixError(MatrixError), } -pub fn sample, const D: usize, #[cfg(feature = "log")] L: Logger>( +pub fn sample( tropical_subgraph_table: &TropicalSubgraphTable, x_space_point: &[T], loop_signature: &[Vec], @@ -31,6 +31,8 @@ pub fn sample, const D: usize, #[cfg(feature = "log")] let num_loops = tropical_subgraph_table.tropical_graph.num_loops; let mut mimic_rng = MimicRng::new(x_space_point); + let const_builder = mimic_rng.zero(); + let permatuhedral_sample = permatuhedral_sampling( tropical_subgraph_table, &mut mimic_rng, @@ -47,27 +49,27 @@ pub fn sample, const D: usize, #[cfg(feature = "log")] } }; - let lambda = Into::::into(inverse_gamma_lr( - tropical_subgraph_table.tropical_graph.dod, - Into::::into(mimic_rng.get_random_number(Some("sample lambda"))), + let lambda = inverse_gamma_lr( + &const_builder.from_f64(tropical_subgraph_table.tropical_graph.dod), + mimic_rng.get_random_number(Some("sample lambda")), 50, - 5.0, - )); + &const_builder.from_f64(5.0), + ); if settings.print_debug_info { #[cfg(not(feature = "log"))] - println!("lambda: {}", lambda); + println!("lambda: {:?}", lambda); #[cfg(feature = "log")] logger.write("momtrop_lambda", &lambda.into()) } - let (edge_masses, edge_shifts): (Vec, Vec>) = edge_data + let (edge_masses, edge_shifts): (Vec, Vec<&Vector>) = edge_data .iter() .map(|(option_mass, edge_shift)| { if let Some(mass) = option_mass { - (*mass, edge_shift) + (mass.clone(), edge_shift) } else { - (T::zero(), edge_shift) + (const_builder.zero(), edge_shift) } }) .unzip(); @@ -84,8 +86,8 @@ pub fn sample, const D: usize, #[cfg(feature = "log")] ); let loop_momenta = compute_loop_momenta( - v_polynomial, - lambda, + &v_polynomial, + &lambda, &decomposed_l_matrix.q_transposed_inverse, &q_vectors, &decomposed_l_matrix.inverse, @@ -95,8 +97,8 @@ pub fn sample, const D: usize, #[cfg(feature = "log")] if settings.print_debug_info { #[cfg(not(feature = "log"))] { - println!("v: {}", v_polynomial); - println!("u: {}", decomposed_l_matrix.determinant); + println!("v: {:?}", &v_polynomial); + println!("u: {:?}", &decomposed_l_matrix.determinant); } #[cfg(feature = "log")] { @@ -107,14 +109,14 @@ pub fn sample, const D: usize, #[cfg(feature = "log")] let u_trop = permatuhedral_sample.u_trop; let v_trop = permatuhedral_sample.v_trop; - let u = decomposed_l_matrix.determinant; + let u = &decomposed_l_matrix.determinant; let v = v_polynomial; - let jacobian = (u_trop / u).powf(Into::::into( - tropical_subgraph_table.dimension as f64 / 2.0, - )) * (v_trop / v) - .powf(Into::::into(tropical_subgraph_table.tropical_graph.dod)) - * Into::::into(tropical_subgraph_table.cached_factor); + let jacobian = (u_trop.ref_div(u)) + .powf(&const_builder.from_f64(tropical_subgraph_table.dimension as f64 / 2.0)) + * (v_trop.ref_div(&v)) + .powf(&const_builder.from_f64(tropical_subgraph_table.tropical_graph.dod)) + * const_builder.from_f64(tropical_subgraph_table.cached_factor); let metadata = if settings.return_metadata { Some(Metadata { @@ -122,7 +124,7 @@ pub fn sample, const D: usize, #[cfg(feature = "log")] q_vectors, lambda, shift: compute_only_shift(&decomposed_l_matrix.inverse, &u_vectors), - decompoisiton_result: decomposed_l_matrix, + decompoisiton_result: decomposed_l_matrix.clone(), u_vectors, }) } else { @@ -133,14 +135,14 @@ pub fn sample, const D: usize, #[cfg(feature = "log")] loop_momenta, u_trop, v_trop, - u, + u: decomposed_l_matrix.determinant, v, jacobian, metadata, }) } -struct PermatuhedralSamplingResult { +struct PermatuhedralSamplingResult { x: Vec, u_trop: T, v_trop: T, @@ -149,16 +151,16 @@ struct PermatuhedralSamplingResult { /// This function returns the feynman parameters for a given graph and sample point, it also computes u_trop and v_trop. /// A rescaling is performed for numerical stability, with this rescaling u_trop and v_trop always evaluate to 1. #[inline] -fn permatuhedral_sampling( +fn permatuhedral_sampling( tropical_subgraph_table: &TropicalSubgraphTable, rng: &mut MimicRng, settings: &TropicalSamplingSettings, #[cfg(feature = "log")] logger: &L, ) -> PermatuhedralSamplingResult { - let mut kappa = T::one(); - let mut x_vec = vec![T::zero(); tropical_subgraph_table.tropical_graph.topology.len()]; - let mut u_trop = T::one(); - let mut v_trop = T::one(); + let mut kappa = rng.one(); + let mut x_vec = vec![rng.zero(); tropical_subgraph_table.tropical_graph.topology.len()]; + let mut u_trop = rng.one(); + let mut v_trop = rng.one(); let mut graph = tropical_subgraph_table .tropical_graph @@ -177,18 +179,18 @@ fn permatuhedral_sampling( tropical_subgraph_table.sample_edge(rng.get_random_number(Some("sample_edge")), &graph) }; - x_vec[edge] = kappa; + x_vec[edge] = kappa.clone(); if tropical_subgraph_table.table[graph.get_id()].mass_momentum_spanning && !tropical_subgraph_table.table[graph_without_edge.get_id()].mass_momentum_spanning { - v_trop = x_vec[edge]; + v_trop = x_vec[edge].clone(); } if tropical_subgraph_table.table[graph_without_edge.get_id()].loop_number < tropical_subgraph_table.table[graph.get_id()].loop_number { - u_trop *= x_vec[edge]; + u_trop *= &x_vec[edge]; } // Terminate early, so we do not waste a random variable in the final step @@ -198,26 +200,27 @@ fn permatuhedral_sampling( } let xi = rng.get_random_number(Some("sample xi")); - kappa *= xi.powf( - Into::::into(tropical_subgraph_table.table[graph.get_id()].generalized_dod).inv(), + kappa *= &xi.powf( + &xi.from_f64(tropical_subgraph_table.table[graph.get_id()].generalized_dod) + .inv(), ); } - let xi_trop = u_trop * v_trop; + let xi_trop = u_trop.ref_mul(&v_trop); // perform rescaling for numerical stability - let target = u_trop.powf(Into::::into( - -(tropical_subgraph_table.dimension as f64 / 2.0), - )) * (u_trop / xi_trop) - .powf(Into::::into(tropical_subgraph_table.tropical_graph.dod)); + let target = u_trop.powf(&xi_trop.from_f64(-(tropical_subgraph_table.dimension as f64 / 2.0))) + * (u_trop.ref_div(&xi_trop)) + .powf(&xi_trop.from_f64(tropical_subgraph_table.tropical_graph.dod)); let loop_number = tropical_subgraph_table.table.last().unwrap().loop_number; let scaling = target.powf( - Into::::into( - tropical_subgraph_table.dimension as f64 / 2.0 * loop_number as f64 - + tropical_subgraph_table.tropical_graph.dod, - ) - .inv(), + &xi_trop + .from_f64( + tropical_subgraph_table.dimension as f64 / 2.0 * loop_number as f64 + + tropical_subgraph_table.tropical_graph.dod, + ) + .inv(), ); if settings.print_debug_info { @@ -230,14 +233,14 @@ fn permatuhedral_sampling( ); } - x_vec.iter_mut().for_each(|x| *x *= scaling); + x_vec.iter_mut().for_each(|x| *x *= &scaling); if settings.print_debug_info { #[cfg(not(feature = "log"))] { - println!("sampled feynman parameters: {:?}", x_vec); - println!("u_trop before rescaling: {:+16e}", u_trop); - println!("v_trop before rescaling: {:+16e}", v_trop); + println!("sampled feynman parameters: {:?}", &x_vec); + println!("u_trop before rescaling: {:?}", &u_trop); + println!("v_trop before rescaling: {:?}", &v_trop); } #[cfg(feature = "log")] { @@ -250,8 +253,8 @@ fn permatuhedral_sampling( } } - u_trop = T::one(); - v_trop = T::one(); + u_trop = u_trop.one(); + v_trop = v_trop.one(); PermatuhedralSamplingResult { x: x_vec, @@ -262,22 +265,25 @@ fn permatuhedral_sampling( /// Compute the L x L matrix from the feynman parameters and the signature matrix #[inline] -fn compute_l_matrix(x_vec: &[T], signature_matrix: &[Vec]) -> SquareMatrix { +fn compute_l_matrix( + x_vec: &[T], + signature_matrix: &[Vec], +) -> SquareMatrix { let num_edges = signature_matrix.len(); let num_loops = signature_matrix[0].len(); - let mut temp_l_matrix = SquareMatrix::new_zeros(num_loops); + let mut temp_l_matrix = SquareMatrix::new_zeros_from_num(&x_vec[0], num_loops); for i in 0..num_loops { for j in i..num_loops { for e in 0..num_edges { - let add = x_vec[e] - * Into::::into((signature_matrix[e][i] * signature_matrix[e][j]) as f64); + let add = x_vec[e].from_isize(signature_matrix[e][i] * signature_matrix[e][j]) + * &x_vec[e]; if i == j { - temp_l_matrix[(i, j)] += add; + temp_l_matrix[(i, j)] += &add; } else { - temp_l_matrix[(i, j)] += add; - temp_l_matrix[(j, i)] += add; + temp_l_matrix[(i, j)] += &add; + temp_l_matrix[(j, i)] += &add; } } } @@ -288,7 +294,7 @@ fn compute_l_matrix(x_vec: &[T], signature_matrix: &[Vec]) /// Sample Gaussian distributed vectors, using the Box-Muller transform #[inline] -fn sample_q_vectors( +fn sample_q_vectors( rng: &mut MimicRng, dimension: usize, num_loops: usize, @@ -296,6 +302,8 @@ fn sample_q_vectors( let token = Some("box muller"); let num_variables = dimension * num_loops; + let builder = rng.zero(); + let num_uniform_variables = num_variables + num_variables % 2; let mut gaussians = (0..num_uniform_variables / 2).flat_map(|_| { let (box_muller_1, box_muller_2) = @@ -307,7 +315,7 @@ fn sample_q_vectors( let mut res = Vec::with_capacity(num_loops); for _ in 0..num_loops { - let mut vec = Vector::::new(); + let mut vec = Vector::::new_from_num(&builder); for i in 0..D { vec[i] = gaussians.next().unwrap_or_else(|| unreachable!()); } @@ -319,46 +327,53 @@ fn sample_q_vectors( /// Compute the vectors u, according to the formula in the notes #[inline] -fn compute_u_vectors( +fn compute_u_vectors( x_vec: &[T], signature_marix: &[Vec], - edge_shifts: &[Vector], + edge_shifts: &[&Vector], ) -> Vec> { let num_loops = signature_marix[0].len(); let num_edges = signature_marix.len(); + let const_builder = &x_vec[0]; (0..num_loops) .map(|l| { - (0..num_edges).fold(Vector::new(), |acc: Vector, e| { - &acc + &(&edge_shifts[e] - * (x_vec[e] * Into::::into(signature_marix[e][l] as f64))) - }) + (0..num_edges).fold( + Vector::new_from_num(const_builder), + |acc: Vector, e| { + &acc + &(edge_shifts[e] + * (const_builder.from_isize(signature_marix[e][l]) * &x_vec[e])) + }, + ) }) .collect_vec() } /// Compute the polynomial v, according to the formula in the notes #[inline] -fn compute_v_polynomial( +fn compute_v_polynomial( x_vec: &[T], u_vectors: &[Vector], inverse_l: &SquareMatrix, - edge_shifts: &[Vector], + edge_shifts: &[&Vector], edge_masses: &[T], ) -> T { let num_loops = inverse_l.get_dim(); + let const_builder = &x_vec[0]; let mut res = izip!(x_vec, edge_masses, edge_shifts) - .map(|(&x_e, &mass, shift)| x_e * (mass * mass + shift.squared())) - .sum::(); + .map(|(x_e, mass, shift)| (mass.ref_mul(mass) + shift.squared()) * x_e) + .fold(const_builder.zero(), |acc, x| acc + x); for l in 0..num_loops { - res -= u_vectors[l].squared() * inverse_l[(l, l)]; + res -= &(u_vectors[l].squared() * &inverse_l[(l, l)]); } for i in 0..num_loops { for j in i + 1..num_loops { - res -= Into::::into(2.) * u_vectors[i].dot(&u_vectors[j]) * inverse_l[(i, j)]; + res -= &(const_builder.from_isize(2) + * u_vectors[i].dot(&u_vectors[j]) + * &inverse_l[(i, j)]); } } @@ -367,24 +382,24 @@ fn compute_v_polynomial( /// Compute the loop momenta, according to the formula in the notes #[inline] -fn compute_loop_momenta( - v: T, - lambda: T, +fn compute_loop_momenta( + v: &T, + lambda: &T, q_t_inverse: &SquareMatrix, q_vectors: &[Vector], l_inverse: &SquareMatrix, u_vectors: &[Vector], ) -> Vec> { let num_loops = q_t_inverse.get_dim(); - let prefactor = (v / lambda * Into::::into(0.5)).sqrt(); + let prefactor = (v.ref_div(lambda) / lambda.from_isize(2)).sqrt(); (0..num_loops) .map(|l| { q_vectors.iter().zip(u_vectors.iter()).enumerate().fold( - Vector::new(), + q_vectors[0].new(), |acc, (l_prime, (q, u))| { - let q_part: Vector = q * (prefactor * q_t_inverse[(l, l_prime)]); - let u_part: Vector = u * l_inverse[(l, l_prime)]; + let q_part: Vector = q * (prefactor.ref_mul(&q_t_inverse[(l, l_prime)])); + let u_part: Vector = u * &l_inverse[(l, l_prime)]; &(&acc + &q_part) - &u_part }, @@ -393,20 +408,21 @@ fn compute_loop_momenta( .collect_vec() } -fn compute_only_shift( +fn compute_only_shift( l_inverse: &SquareMatrix, u_vectors: &[Vector], ) -> Vec> { let num_loops = l_inverse.get_dim(); + (0..num_loops) .map(|l| { - u_vectors - .iter() - .enumerate() - .fold(Vector::new(), |acc, (l_prime, u)| { - let u_part: Vector = u * l_inverse[(l, l_prime)]; + u_vectors.iter().enumerate().fold( + Vector::new_from_num(&l_inverse.zero()), + |acc, (l_prime, u)| { + let u_part: Vector = u * &l_inverse[(l, l_prime)]; &acc + &u_part - }) + }, + ) }) .collect_vec() } diff --git a/src/vector.rs b/src/vector.rs index 7292673..c966d32 100644 --- a/src/vector.rs +++ b/src/vector.rs @@ -1,13 +1,15 @@ -use crate::float::FloatLike; -use f128::f128; -use std::ops::{Add, AddAssign, Index, IndexMut, Mul, Sub}; +use crate::float::MomTropFloat; +use std::{ + array, + ops::{Add, AddAssign, Index, IndexMut, Mul, Sub}, +}; #[derive(Clone, Copy, Debug)] -pub struct Vector { +pub struct Vector { elements: [T; D], } -impl Vector { +impl Vector { #[inline] pub fn from_array(elements: [T; D]) -> Self { Self { elements } @@ -22,39 +24,47 @@ impl Vector { } } + #[inline] + pub fn zero(&self) -> T { + self.elements[0].zero() + } + #[inline] pub fn from_slice(elements: &[T; D]) -> Self { Self { - elements: *elements, + elements: elements.clone(), } } #[inline] - pub fn new() -> Self { + pub fn new(&self) -> Self { Self { - elements: [T::zero(); D], + elements: self.elements.each_ref().map(|value| value.zero()), } } #[inline] - pub fn squared(&self) -> T { - let mut res = T::zero(); - for elem in self.elements { - res += elem * elem + pub fn new_from_num(builder: &T) -> Self { + Self { + elements: array::from_fn(|_| builder.zero()), } + } - res + #[inline] + pub fn squared(&self) -> T { + self.elements + .iter() + .fold(self.elements[0].zero(), |acc, x| acc + x.ref_mul(x)) } #[inline] pub fn dot(&self, rhs: &Self) -> T { - let mut res = T::zero(); - - for (&lhs_elem, &rhs_elem) in self.elements.iter().zip(rhs.elements.iter()) { - res += lhs_elem * rhs_elem - } - - res + self.elements + .iter() + .zip(rhs.elements.iter()) + .fold(self.elements[0].zero(), |acc, (left, right)| { + acc + left.ref_mul(right) + }) } // We are not going to do zero-dimensional qft! @@ -66,11 +76,11 @@ impl Vector { #[inline] pub fn get_elements(&self) -> [T; D] { - self.elements + self.elements.clone() } } -impl Index for Vector { +impl Index for Vector { type Output = T; fn index(&self, index: usize) -> &Self::Output { @@ -78,91 +88,88 @@ impl Index for Vector { } } -impl IndexMut for Vector { +impl IndexMut for Vector { fn index_mut(&mut self, index: usize) -> &mut Self::Output { &mut self.elements[index] } } -impl Vector { - pub fn upcast(&self) -> Vector { - let mut elements = [f128::new(0.0); D]; - for (new_element, old_element) in elements.iter_mut().zip(self.elements.iter()) { - *new_element = f128::new(*old_element); - } +impl Mul for &Vector { + type Output = Vector; - Vector { elements } + fn mul(self, rhs: T) -> Self::Output { + Self::Output { + elements: self.elements.each_ref().map(|elem| elem.ref_mul(&rhs)), + } } } -impl Mul for &Vector { +impl Mul<&T> for &Vector { type Output = Vector; - fn mul(self, rhs: T) -> Self::Output { - let mut res = *self; - - for i in 0..D { - res[i] *= rhs; + fn mul(self, rhs: &T) -> Self::Output { + Self::Output { + elements: self.elements.each_ref().map(|elem| elem.ref_mul(rhs)), } - - res } } -impl Add<&Vector> for &Vector { +impl Add<&Vector> for &Vector { type Output = Vector; #[inline] fn add(self, rhs: &Vector) -> Self::Output { - let mut res = Self::Output::new(); - - for i in 0..D { - res[i] += self[i] + rhs[i]; + //Self::Output { + // elements: self + // .elements + // .iter() + // .zip(&rhs.elements) + // .map(|(left, right)| left + right) + // .collect_vec() + // .try_into() + // .unwrap_or_else(|_| unreachable!()), + //} + + Self::Output { + elements: array::from_fn(|i| self[i].ref_add(&rhs[i])), } - - res } } -impl Sub<&Vector> for &Vector { +impl Sub<&Vector> for &Vector { type Output = Vector; #[inline] fn sub(self, rhs: &Vector) -> Self::Output { - let mut res = Self::Output::new(); - - for i in 0..D { - res[i] += self[i] - rhs[i]; + //Self::Output { + // elements: self + // .elements + // .iter() + // .zip(&rhs.elements) + // .map(|(left, right)| left - right) + // .collect_vec() + // .try_into() + // .unwrap_or_else(|_| unreachable!()), + //} + + Self::Output { + elements: array::from_fn(|i| self[i].ref_sub(&rhs[i])), } - - res } } -impl AddAssign for Vector { +impl AddAssign for Vector { #[inline] fn add_assign(&mut self, rhs: Self) { for i in 0..D { - self[i] += rhs[i]; - } - } -} - -impl Vector { - pub fn downcast(&self) -> Vector { - let mut new_elements: [f64; D] = [0.0; D]; - for (new_element, out_element) in new_elements.iter_mut().zip(self.elements) { - *new_element = out_element.into(); - } - - Vector { - elements: new_elements, + self[i] += &rhs[i]; } } } #[cfg(test)] mod tests { + use super::Vector; #[test] @@ -187,7 +194,7 @@ mod tests { #[test] fn test_new() { - let vector: Vector = Vector::new(); + let vector = Vector::from_array([1., 2., 3.]).new(); assert_eq!(vector.elements[0], 0.0); assert_eq!(vector.elements[1], 0.0); @@ -209,15 +216,6 @@ mod tests { assert_eq!(dot, 0.0); } - #[test] - fn test_upcast() { - let vector = Vector::from_array([2.0, 3.0]); - let vector_f128 = vector.upcast(); - - assert_eq!(vector_f128.elements[0], f128::f128::new(2.0)); - assert_eq!(vector_f128.elements[1], f128::f128::new(3.0)); - } - #[test] fn test_mul() { let vector = Vector::from_array([2.0, 4.0]); diff --git a/tests/test_prec_eq.rs b/tests/test_prec_eq.rs index 13bd43a..e64d016 100644 --- a/tests/test_prec_eq.rs +++ b/tests/test_prec_eq.rs @@ -1,3 +1,4 @@ +use momtrop::assert_approx_eq; use momtrop::{ vector::Vector, Edge, Graph, SampleGenerator, TropicalSampleResult, TropicalSamplingSettings, }; @@ -37,12 +38,12 @@ fn test_prec_eq() { let p2 = Vector::from_array([6.0, 7.0, 8.0]); let edge_data = vec![ - (Option::::None, Vector::new()), + (Option::::None, p1.new()), (None, p1), (None, (&p1 + &p2)), ]; - let x_space_point = vec![0.1; sampler.get_dimension()]; + let x_space_point = vec![(0.1); sampler.get_dimension()]; let settings = TropicalSamplingSettings::default(); #[cfg(feature = "log")] @@ -57,20 +58,21 @@ fn test_prec_eq() { &logger, ) .unwrap(); - let sample_f128 = sampler - .generate_sample_f128_from_x_space_point( - &x_space_point, - edge_data.clone(), - &settings, - #[cfg(feature = "log")] - &logger, - ) - .unwrap() - .downcast(); + //let sample_f128 = sampler + // .generate_sample_f128_from_x_space_point( + // &x_space_point, + // edge_data.clone(), + // &settings, + // #[cfg(feature = "log")] + // &logger, + // ) + // .unwrap() + // .downcast(); - assert_approx_eq_sample(sample, sample_f128, 1.0e-15); + //assert_approx_eq_sample(sample, sample_f128, 1.0e-15); } +#[cfg(test)] fn assert_approx_eq_sample( sample_1: TropicalSampleResult, sample_2: TropicalSampleResult, @@ -84,20 +86,14 @@ fn assert_approx_eq_sample( let elements_1 = loop_mom_1.get_elements(); let elements_2 = loop_mom_2.get_elements(); - for (&el_1, &el_2) in elements_1.iter().zip(elements_2.iter()) { - assert_approx_eq(el_1, el_2, tolerance); + for (el_1, el_2) in elements_1.iter().zip(elements_2.iter()) { + assert_approx_eq(el_1, el_2, &tolerance); } } - assert_approx_eq(sample_1.u_trop, sample_2.u_trop, tolerance); - assert_approx_eq(sample_1.v_trop, sample_2.v_trop, tolerance); - assert_approx_eq(sample_1.u, sample_2.u, tolerance); - assert_approx_eq(sample_1.v, sample_2.v, tolerance); - assert_approx_eq(sample_1.jacobian, sample_2.jacobian, tolerance); -} - -fn assert_approx_eq(x: f64, y: f64, tolerance: f64) { - let avg = (x + y) / 2.; - let norm_err = ((x - y) / avg).abs(); - assert!(norm_err < tolerance); + assert_approx_eq(&sample_1.u_trop, &sample_2.u_trop, &tolerance); + assert_approx_eq(&sample_1.v_trop, &sample_2.v_trop, &tolerance); + assert_approx_eq(&sample_1.u, &sample_2.u, &tolerance); + assert_approx_eq(&sample_1.v, &sample_2.v, &tolerance); + assert_approx_eq(&sample_1.jacobian, &sample_2.jacobian, &tolerance); } diff --git a/tests/triangle.rs b/tests/triangle.rs index 88b9544..378e498 100644 --- a/tests/triangle.rs +++ b/tests/triangle.rs @@ -1,4 +1,4 @@ -use momtrop::{vector::Vector, Edge, Graph, TropicalSamplingSettings}; +use momtrop::{float::MomTropFloat, vector::Vector, Edge, Graph, TropicalSamplingSettings}; use rand::SeedableRng; /// integrate a massless triangle with LTD and tropicalsampling @@ -44,13 +44,18 @@ fn integrate_massless_triangle() { let mut min_pol_ratio: f64 = 1.0; let n_samples = 1000000; - - let edge_data = vec![(None, Vector::new()), (None, p1), (None, (&p1 + &p2))]; - let p10 = 1.0; let p20 = 1.0; - let settings = TropicalSamplingSettings::default(); + let edge_data = vec![ + (None, Vector::new_from_num(&p10)), + (None, p1), + (None, (&p1 + &p2)), + ]; + + let settings = TropicalSamplingSettings { + ..Default::default() + }; #[cfg(feature = "log")] let logger = momtrop::log::DummyLogger {}; @@ -78,31 +83,32 @@ fn integrate_massless_triangle() { max_pol_ratio = max_pol_ratio.max(polynomial_ratio); min_pol_ratio = min_pol_ratio.min(polynomial_ratio); - let pi_factor = (2f64 * std::f64::consts::PI).powi(3); + let pi_factor = (p10.from_isize(2) * p10.PI()).powf(p10.from_isize(3)); let prefactor = energy_prefactor * sample.jacobian / pi_factor; - let term1 = ((energy_0 + energy_1 + p10) * (energy_2 + energy_0 + p10 + p20)).recip(); - let term2 = ((energy_2 + energy_0 - p10 - p20) * (energy_1 + energy_2 - p20)).recip(); - let term3 = ((energy_0 + energy_1 + p10) * (energy_1 + energy_2 - p20)).recip(); + let term1 = ((energy_0 + energy_1 + p10) * (energy_2 + energy_0 + p10 + p20)).inv(); + let term2 = ((energy_2 + energy_0 - p10 - p20) * (energy_1 + energy_2 - p20)).inv(); + let term3 = ((energy_0 + energy_1 + p10) * (energy_1 + energy_2 - p20)).inv(); - let term4 = ((energy_0 + energy_1 - p10) * (energy_2 + energy_0 - p10 - p20)).recip(); - let term5 = ((energy_2 + energy_0 + p10 + p20) * (energy_1 + energy_2 + p20)).recip(); - let term6 = ((energy_0 + energy_1 - p10) * (energy_1 + energy_2 + p20)).recip(); + let term4 = ((energy_0 + energy_1 - p10) * (energy_2 + energy_0 - p10 - p20)).inv(); + let term5 = ((energy_2 + energy_0 + p10 + p20) * (energy_1 + energy_2 + p20)).inv(); + let term6 = ((energy_0 + energy_1 - p10) * (energy_1 + energy_2 + p20)).inv(); sum += (term1 + term2 + term3 + term4 + term5 + term6) * prefactor; } - let avg = sum / n_samples as f64; + let avg = sum / p10.from_isize(n_samples as isize); // theoretical bound assert!( min_pol_ratio - >= (1f64 / 3f64).powf(sampler.get_dimension() as f64 / 2. - sampler.get_dod()) - * (1f64 / (p1.squared() + p2.squared() + (&p1 + &p2).squared())) - .powf(sampler.get_dod()) + >= (p10.one() / p10.from_isize(3)) + .powf(sampler.get_dimension() as f64 / 2. - sampler.get_dod()) + * (p10.one() / (p1.squared() + p2.squared() + (&p1 + &p2).squared())) + .powf(p10.from_f64(sampler.get_dod())) ); - assert!(max_pol_ratio <= (1f64 / p1.squared()).powf(sampler.get_dod())); + assert!(max_pol_ratio <= (p10.one() / p1.squared()).powf(sampler.get_dod())); // this is the exact value with this seed, needs a more robust test assert_eq!(9.758362839019336e-5, avg); From d4061a9ab8175b324edc914fabae05cfbe1af48b Mon Sep 17 00:00:00 2001 From: Mathijs Fraaije Date: Wed, 27 Nov 2024 11:45:50 +0100 Subject: [PATCH 2/6] fix log feature --- src/gamma.rs | 39 ++++++++++++++++++++++++++++++--------- src/sampling.rs | 14 +++++++------- 2 files changed, 37 insertions(+), 16 deletions(-) diff --git a/src/gamma.rs b/src/gamma.rs index 2cda048..040e001 100644 --- a/src/gamma.rs +++ b/src/gamma.rs @@ -1,7 +1,24 @@ use statrs::function::gamma::{gamma, gamma_lr, gamma_ur}; +use crate::float::MomTropFloat; + +#[inline] +pub fn inverse_gamma_lr( + a: &T, + p: &T, + max_n_iter: usize, + epsilon_tolerance: &T, +) -> T { + a.from_f64(inverse_gamma_lr_impl( + a.to_f64(), + p.to_f64(), + max_n_iter, + epsilon_tolerance.to_f64(), + )) +} + #[inline] -pub fn inverse_gamma_lr(a: f64, p: f64, max_n_iter: usize, epsilon_tolerance: f64) -> f64 { +pub fn inverse_gamma_lr_impl(a: f64, p: f64, max_n_iter: usize, epsilon_tolerance: f64) -> f64 { // this algorithm is taken from https://dl.acm.org/doi/pdf/10.1145/22721.23109 // get an estimate for x0 to start newton iterations. @@ -189,8 +206,9 @@ mod tests { let omega = 1.; for p in [0.1, 0.3, 0.5, 0.7, 0.9] { - let inverse_lower_gamma = inverse_gamma_lr(omega, p, NITER_FOR_TEST, EPSILON_TOLERANCE); - assert_approx_eq(gamma_lr(omega, inverse_lower_gamma), p, TOLERANCE); + let inverse_lower_gamma = + inverse_gamma_lr(&omega, &p, NITER_FOR_TEST, &EPSILON_TOLERANCE); + assert_approx_eq(&(gamma_lr(omega, inverse_lower_gamma)), &p, &TOLERANCE); } } @@ -198,8 +216,9 @@ mod tests { fn test_half() { let omega = 0.5; for p in [0.1, 0.3, 0.5, 0.7, 0.9] { - let inverse_lower_gamma = inverse_gamma_lr(omega, p, NITER_FOR_TEST, EPSILON_TOLERANCE); - assert_approx_eq(gamma_lr(omega, inverse_lower_gamma), p, TOLERANCE); + let inverse_lower_gamma = + inverse_gamma_lr(&omega, &p, NITER_FOR_TEST, &EPSILON_TOLERANCE); + assert_approx_eq(&(gamma_lr(omega, inverse_lower_gamma)), &p, &TOLERANCE); } } @@ -207,8 +226,9 @@ mod tests { fn test_2() { let omega = 2.0; for p in [0.1, 0.3, 0.5, 0.7, 0.9] { - let inverse_lower_gamma = inverse_gamma_lr(omega, p, NITER_FOR_TEST, EPSILON_TOLERANCE); - assert_approx_eq(gamma_lr(omega, inverse_lower_gamma), p, TOLERANCE); + let inverse_lower_gamma = + inverse_gamma_lr(&omega, &p, NITER_FOR_TEST, &EPSILON_TOLERANCE); + assert_approx_eq(&(gamma_lr(omega, inverse_lower_gamma)), &p, &TOLERANCE); } } @@ -216,8 +236,9 @@ mod tests { fn test_10() { let omega = 10.0; for p in [0.1, 0.3, 0.5, 0.7, 0.9] { - let inverse_lower_gamma = inverse_gamma_lr(omega, p, NITER_FOR_TEST, EPSILON_TOLERANCE); - assert_approx_eq(gamma_lr(omega, inverse_lower_gamma), p, TOLERANCE); + let inverse_lower_gamma = + inverse_gamma_lr(&omega, &p, NITER_FOR_TEST, &EPSILON_TOLERANCE); + assert_approx_eq(&(gamma_lr(omega, inverse_lower_gamma)), &p, &TOLERANCE); } } } diff --git a/src/sampling.rs b/src/sampling.rs index fd6fb76..d80da67 100644 --- a/src/sampling.rs +++ b/src/sampling.rs @@ -60,7 +60,7 @@ pub fn sample, Vec<&Vector>) = edge_data @@ -102,8 +102,8 @@ pub fn sample( #[cfg(feature = "log")] logger.write( "momtrop_feynman_parameter_no_rescaling", - &x_vec.iter().map(|x| Into::::into(*x)).collect_vec(), + &x_vec.iter().map(|x| x.to_f64()).collect_vec(), ); } @@ -246,10 +246,10 @@ fn permatuhedral_sampling( { logger.write( "momtrop_feynman_parameter", - &x_vec.iter().map(|x| Into::::into(*x)).collect_vec(), + &x_vec.iter().map(|x| x.to_f64()).collect_vec(), ); - logger.write("momtrop_u_trop_no_rescaling", &u_trop.into()); - logger.write("momtrop_v_trop_no_rescaling", &v_trop.into()); + logger.write("momtrop_u_trop_no_rescaling", &u_trop.to_f64()); + logger.write("momtrop_v_trop_no_rescaling", &v_trop.to_f64()); } } From e91828b93fc777613724d5e57f6a0a1a4d450027 Mon Sep 17 00:00:00 2001 From: Mathijs Fraaije Date: Wed, 4 Dec 2024 11:35:25 +0100 Subject: [PATCH 3/6] make lambda sampling return error if it generates a NaN --- src/gamma.rs | 23 ++++++++++++++++------- src/sampling.rs | 5 ++++- 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/src/gamma.rs b/src/gamma.rs index 040e001..030aa2b 100644 --- a/src/gamma.rs +++ b/src/gamma.rs @@ -2,19 +2,28 @@ use statrs::function::gamma::{gamma, gamma_lr, gamma_ur}; use crate::float::MomTropFloat; +#[derive(Debug, Clone, Copy)] +pub struct GammaError {} + #[inline] pub fn inverse_gamma_lr( a: &T, p: &T, max_n_iter: usize, epsilon_tolerance: &T, -) -> T { - a.from_f64(inverse_gamma_lr_impl( +) -> Result { + let res = inverse_gamma_lr_impl( a.to_f64(), p.to_f64(), max_n_iter, epsilon_tolerance.to_f64(), - )) + ); + + if res.is_nan() { + Err(GammaError {}) + } else { + Ok(a.from_f64(res)) + } } #[inline] @@ -207,7 +216,7 @@ mod tests { for p in [0.1, 0.3, 0.5, 0.7, 0.9] { let inverse_lower_gamma = - inverse_gamma_lr(&omega, &p, NITER_FOR_TEST, &EPSILON_TOLERANCE); + inverse_gamma_lr(&omega, &p, NITER_FOR_TEST, &EPSILON_TOLERANCE).unwrap(); assert_approx_eq(&(gamma_lr(omega, inverse_lower_gamma)), &p, &TOLERANCE); } } @@ -217,7 +226,7 @@ mod tests { let omega = 0.5; for p in [0.1, 0.3, 0.5, 0.7, 0.9] { let inverse_lower_gamma = - inverse_gamma_lr(&omega, &p, NITER_FOR_TEST, &EPSILON_TOLERANCE); + inverse_gamma_lr(&omega, &p, NITER_FOR_TEST, &EPSILON_TOLERANCE).unwrap(); assert_approx_eq(&(gamma_lr(omega, inverse_lower_gamma)), &p, &TOLERANCE); } } @@ -227,7 +236,7 @@ mod tests { let omega = 2.0; for p in [0.1, 0.3, 0.5, 0.7, 0.9] { let inverse_lower_gamma = - inverse_gamma_lr(&omega, &p, NITER_FOR_TEST, &EPSILON_TOLERANCE); + inverse_gamma_lr(&omega, &p, NITER_FOR_TEST, &EPSILON_TOLERANCE).unwrap(); assert_approx_eq(&(gamma_lr(omega, inverse_lower_gamma)), &p, &TOLERANCE); } } @@ -237,7 +246,7 @@ mod tests { let omega = 10.0; for p in [0.1, 0.3, 0.5, 0.7, 0.9] { let inverse_lower_gamma = - inverse_gamma_lr(&omega, &p, NITER_FOR_TEST, &EPSILON_TOLERANCE); + inverse_gamma_lr(&omega, &p, NITER_FOR_TEST, &EPSILON_TOLERANCE).unwrap(); assert_approx_eq(&(gamma_lr(omega, inverse_lower_gamma)), &p, &TOLERANCE); } } diff --git a/src/sampling.rs b/src/sampling.rs index d80da67..b7b0a27 100644 --- a/src/sampling.rs +++ b/src/sampling.rs @@ -1,6 +1,7 @@ use itertools::{izip, Itertools}; use super::TropicalSubgraphTable; +use crate::gamma::GammaError; #[cfg(feature = "log")] use crate::log::Logger; use crate::matrix::{MatrixError, SquareMatrix}; @@ -18,6 +19,7 @@ fn box_muller(x1: &T, x2: &T) -> (T, T) { #[derive(Debug, Clone, Copy)] pub enum SamplingError { MatrixError(MatrixError), + GammaError(GammaError), } pub fn sample( @@ -54,7 +56,8 @@ pub fn sample Date: Mon, 13 Jan 2025 12:46:21 +0100 Subject: [PATCH 4/6] change dimension argument to const generic --- benches/graph_bench.rs | 2 +- src/lib.rs | 5 ++--- tests/test_prec_eq.rs | 4 ++-- tests/triangle.rs | 2 +- 4 files changed, 6 insertions(+), 7 deletions(-) diff --git a/benches/graph_bench.rs b/benches/graph_bench.rs index 80baab6..a321471 100644 --- a/benches/graph_bench.rs +++ b/benches/graph_bench.rs @@ -34,7 +34,7 @@ fn criterion_benchmark(c: &mut Criterion) { let logger = momtrop::log::DummyLogger {}; let loop_signature = vec![vec![1]; 3]; - let sampler = graph.build_sampler(loop_signature, 3).unwrap(); + let sampler = graph.build_sampler(loop_signature).unwrap(); let mut rng = rand::rngs::StdRng::seed_from_u64(69); let p1 = Vector::from_array([3.0, 4.0, 5.0]); let p2 = Vector::from_array([6.0, 7.0, 8.0]); diff --git a/src/lib.rs b/src/lib.rs index d56656d..dc8e917 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -90,10 +90,9 @@ impl Graph { pub fn build_sampler( self, loop_signature: Vec>, - dimension: usize, ) -> Result, String> { - let tropical_graph = TropicalGraph::from_graph(self, dimension); - let table = TropicalSubgraphTable::generate_from_tropical(&tropical_graph, dimension)?; + let tropical_graph = TropicalGraph::from_graph(self, D); + let table = TropicalSubgraphTable::generate_from_tropical(&tropical_graph, D)?; Ok(SampleGenerator { loop_signature, diff --git a/tests/test_prec_eq.rs b/tests/test_prec_eq.rs index e64d016..4ccf02b 100644 --- a/tests/test_prec_eq.rs +++ b/tests/test_prec_eq.rs @@ -33,7 +33,7 @@ fn test_prec_eq() { }; let loop_signature = vec![vec![1]; 3]; - let sampler: SampleGenerator<3> = graph.build_sampler(loop_signature, 3).unwrap(); + let sampler: SampleGenerator<3> = graph.build_sampler(loop_signature).unwrap(); let p1 = Vector::from_array([3.0, 4.0, 5.0]); let p2 = Vector::from_array([6.0, 7.0, 8.0]); @@ -43,7 +43,7 @@ fn test_prec_eq() { (None, (&p1 + &p2)), ]; - let x_space_point = vec![(0.1); sampler.get_dimension()]; + let x_space_point = vec![0.1; sampler.get_dimension()]; let settings = TropicalSamplingSettings::default(); #[cfg(feature = "log")] diff --git a/tests/triangle.rs b/tests/triangle.rs index 378e498..b7ffdfa 100644 --- a/tests/triangle.rs +++ b/tests/triangle.rs @@ -32,7 +32,7 @@ fn integrate_massless_triangle() { }; let loop_signature = vec![vec![1]; 3]; - let sampler = graph.build_sampler(loop_signature, 3).unwrap(); + let sampler = graph.build_sampler(loop_signature).unwrap(); let p1 = Vector::from_array([3.0, 4.0, 5.0]); let p2 = Vector::from_array([6.0, 7.0, 8.0]); From 2a2b8b72143ea8889a05a5153bf6aef52e15fb7a Mon Sep 17 00:00:00 2001 From: Mathijs Fraaije Date: Mon, 13 Jan 2025 13:15:15 +0100 Subject: [PATCH 5/6] remove upcast_on_failure, the user must handle this themselves since f128 is no longer supplied by default --- src/lib.rs | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index dc8e917..0994229 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -27,16 +27,15 @@ pub const MAX_VERTICES: usize = 256; #[derive(Debug)] pub struct TropicalSamplingSettings { - pub upcast_on_failure: bool, pub matrix_stability_test: Option, pub print_debug_info: bool, pub return_metadata: bool, } +#[allow(clippy::derivable_impls)] impl Default for TropicalSamplingSettings { fn default() -> Self { Self { - upcast_on_failure: true, matrix_stability_test: None, print_debug_info: false, return_metadata: false, From 9a3a0f79f454ce1bd25e28148d3a50d9e4bb3b81 Mon Sep 17 00:00:00 2001 From: Mathijs Fraaije Date: Mon, 13 Jan 2025 13:23:59 +0100 Subject: [PATCH 6/6] remove test for precision equality, as it is no longer relevant --- tests/test_prec_eq.rs | 99 ------------------------------------------- 1 file changed, 99 deletions(-) delete mode 100644 tests/test_prec_eq.rs diff --git a/tests/test_prec_eq.rs b/tests/test_prec_eq.rs deleted file mode 100644 index 4ccf02b..0000000 --- a/tests/test_prec_eq.rs +++ /dev/null @@ -1,99 +0,0 @@ -use momtrop::assert_approx_eq; -use momtrop::{ - vector::Vector, Edge, Graph, SampleGenerator, TropicalSampleResult, TropicalSamplingSettings, -}; - -#[test] -fn test_prec_eq() { - let weight = 2. / 3.; - - let triangle_edges = vec![ - Edge { - vertices: (0, 1), - is_massive: false, - weight, - }, - Edge { - vertices: (1, 2), - is_massive: false, - weight, - }, - Edge { - vertices: (2, 0), - is_massive: false, - weight, - }, - ]; - - let externals = vec![0, 1, 2]; - - let graph = Graph { - edges: triangle_edges, - externals, - }; - - let loop_signature = vec![vec![1]; 3]; - let sampler: SampleGenerator<3> = graph.build_sampler(loop_signature).unwrap(); - let p1 = Vector::from_array([3.0, 4.0, 5.0]); - let p2 = Vector::from_array([6.0, 7.0, 8.0]); - - let edge_data = vec![ - (Option::::None, p1.new()), - (None, p1), - (None, (&p1 + &p2)), - ]; - - let x_space_point = vec![0.1; sampler.get_dimension()]; - let settings = TropicalSamplingSettings::default(); - - #[cfg(feature = "log")] - let logger = momtrop::log::DummyLogger {}; - - let sample = sampler - .generate_sample_from_x_space_point( - &x_space_point, - edge_data.clone(), - &settings, - #[cfg(feature = "log")] - &logger, - ) - .unwrap(); - //let sample_f128 = sampler - // .generate_sample_f128_from_x_space_point( - // &x_space_point, - // edge_data.clone(), - // &settings, - // #[cfg(feature = "log")] - // &logger, - // ) - // .unwrap() - // .downcast(); - - //assert_approx_eq_sample(sample, sample_f128, 1.0e-15); -} - -#[cfg(test)] -fn assert_approx_eq_sample( - sample_1: TropicalSampleResult, - sample_2: TropicalSampleResult, - tolerance: f64, -) { - for (loop_mom_1, loop_mom_2) in sample_1 - .loop_momenta - .iter() - .zip(sample_2.loop_momenta.iter()) - { - let elements_1 = loop_mom_1.get_elements(); - let elements_2 = loop_mom_2.get_elements(); - - for (el_1, el_2) in elements_1.iter().zip(elements_2.iter()) { - assert_approx_eq(el_1, el_2, &tolerance); - } - } - - assert_approx_eq(&sample_1.u_trop, &sample_2.u_trop, &tolerance); - assert_approx_eq(&sample_1.v_trop, &sample_2.v_trop, &tolerance); - assert_approx_eq(&sample_1.u, &sample_2.u, &tolerance); - assert_approx_eq(&sample_1.v, &sample_2.v, &tolerance); - assert_approx_eq(&sample_1.jacobian, &sample_2.jacobian, &tolerance); -}