Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rustify ome.u.s.{as1, as2} #385

Merged
merged 48 commits into from
Oct 14, 2024
Merged
Show file tree
Hide file tree
Changes from 47 commits
Commits
Show all changes
48 commits
Select commit Hold shift + click to select a range
56a77e4
starting ome
tgiani May 31, 2024
95aec14
as1. Need some cleaning up
tgiani Jun 19, 2024
91aaf67
spacelike.rs
tgiani Jun 19, 2024
6eacd3b
attempt to implement python hook
tgiani Jun 21, 2024
b20732d
refactoring and adding Lsv for ome
tgiani Jul 1, 2024
637a1c3
Improve quad_ker module
felixhekhorn Jul 17, 2024
56e9062
rust: Move ome changes to patch
felixhekhorn Jul 17, 2024
12bb7fe
rust: Adjust some docs
felixhekhorn Jul 17, 2024
31f5d88
rust: Fix py callback signature
felixhekhorn Jul 17, 2024
865da0c
rust: Unify unravel of tensors
felixhekhorn Jul 17, 2024
7f055b2
Merge branch 'master' into ome_rust
felixhekhorn Jul 18, 2024
13bafbb
rust: Add several fixes
felixhekhorn Jul 18, 2024
0028c9d
rust: Fix doc abbrev
felixhekhorn Jul 18, 2024
880d7aa
rust: Update docs
felixhekhorn Jul 18, 2024
5f4bde6
rust: Set is_ome flag
felixhekhorn Jul 19, 2024
6323f1b
start as2
tgiani Jul 30, 2024
837187b
Merge branch 'master' into ome_rust
felixhekhorn Aug 9, 2024
7c1a1bd
Fix rust patches
felixhekhorn Aug 9, 2024
2b623d8
Merge branch 'master' into ome_rust
felixhekhorn Aug 22, 2024
4b1755b
rust: Fix a clippy warning
felixhekhorn Aug 22, 2024
56ed86c
Merge branch 'master' into ome_rust
felixhekhorn Aug 22, 2024
164c02b
rust: Protect cache properties
felixhekhorn Aug 22, 2024
33574b9
adding Sm2e and Sm2o implementation
tgiani Sep 4, 2024
e79f008
Merge branch 'ome_rust' of github.com:N3PDF/eko into ome_rust
tgiani Sep 4, 2024
e57e891
adding Sm1e, Sm1o, Sm21e, Sm21o implementation
tgiani Sep 4, 2024
70cdc7c
adding Sm3e and Sm3o
tgiani Sep 5, 2024
678efa9
working on as2
tgiani Sep 5, 2024
7ea6da2
Agq
tgiani Sep 5, 2024
429206a
more work on as2
tgiani Sep 9, 2024
f3b1219
rust: Remove leading underscore for used vars
felixhekhorn Sep 10, 2024
3f27e6b
rust: Add NNLO OME refs
felixhekhorn Sep 10, 2024
17a73d5
rust: Update NLO OME tests
felixhekhorn Sep 10, 2024
5239d8b
rust: Use cmplx eq in tests
felixhekhorn Sep 10, 2024
9b498b7
rust: Simplify assert_approx_eq_cmplx
felixhekhorn Sep 10, 2024
8f09ea0
rust: Activate NNLO OME
felixhekhorn Sep 10, 2024
2d9b7ca
start unit tests on as2
tgiani Sep 11, 2024
fc01791
Merge branch 'ome_rust' of github.com:N3PDF/eko into ome_rust
tgiani Sep 11, 2024
b7932aa
fix bug in Ahqps
tgiani Sep 11, 2024
b1beb83
Merge branch 'master' into ome_rust
felixhekhorn Sep 16, 2024
7d5dac7
Fix wrong merge
felixhekhorn Sep 16, 2024
ddc57bf
Fix OME patch
felixhekhorn Sep 16, 2024
33a1429
tests on as2
tgiani Oct 3, 2024
622f113
more work on as2 tests
tgiani Oct 3, 2024
457e52c
Update crates/ekore/src/operator_matrix_elements/unpolarized/spacelik…
felixhekhorn Oct 11, 2024
d1b6c53
rust: Add rel to cmplx eq
felixhekhorn Oct 11, 2024
7396c67
rust: Unify cmplx! macro call
felixhekhorn Oct 11, 2024
54d86d4
rust: Cast to float with .
felixhekhorn Oct 11, 2024
8896dea
rust: Unify Sm21 in the back
felixhekhorn Oct 14, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions crates/dekoder/src/eko.rs
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ impl EKO {
pub fn write(&self, dst: PathBuf) -> Result<()> {
self.assert_working_dir()?;
// create writer
let dst_file = File::create(&dst)?;
let dst_file = File::create(dst)?;
let dst_file = BufWriter::with_capacity(TAR_WRITER_CAPACITY, dst_file);
let mut ar = tar::Builder::new(dst_file);
// do it!
Expand All @@ -101,7 +101,7 @@ impl EKO {

/// Extract tar file from `src` to `dst`.
pub fn extract(src: PathBuf, dst: PathBuf) -> Result<Self> {
let mut ar = tar::Archive::new(File::open(&src)?);
let mut ar = tar::Archive::new(File::open(src)?);
ar.unpack(&dst)?;
Self::load_opened(dst)
}
Expand Down
97 changes: 75 additions & 22 deletions crates/eko/src/lib.rs
Original file line number Diff line number Diff line change
@@ -1,11 +1,35 @@
//! Interface to the eko Python package.

use ekore::harmonics::cache::Cache;
use num::Complex;
use std::ffi::c_void;

pub mod bib;
pub mod mellin;

/// Wrapper to pass arguments back to Python
struct RawCmplx {
re: Vec<f64>,
im: Vec<f64>,
}

/// Map tensors to c-ordered list
fn unravel<const DIM: usize>(res: Vec<[[Complex<f64>; DIM]; DIM]>, order_qcd: usize) -> RawCmplx {
let mut target = RawCmplx {
re: Vec::<f64>::new(),
im: Vec::<f64>::new(),
};
for obj in res.iter().take(order_qcd) {
for col in obj.iter().take(DIM) {
for el in col.iter().take(DIM) {
target.re.push(el.re);
target.im.push(el.im);
}
}
}
target
}

/// QCD intergration kernel inside quad.
///
/// # Safety
Expand All @@ -14,44 +38,66 @@ pub mod mellin;
pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 {
let args = *(rargs as *mut QuadQCDargs);
let is_singlet = (100 == args.mode0) || (21 == args.mode0) || (90 == args.mode0);
// prepare gamma
// prepare Mellin stuff
let path = mellin::TalbotPath::new(u, args.logx, is_singlet);
let jac = path.jac() * path.prefactor();
let mut c = Cache::new(path.n());
let mut re = Vec::<f64>::new();
let mut im = Vec::<f64>::new();
if is_singlet {
let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd(
let mut raw = RawCmplx {
re: Vec::<f64>::new(),
im: Vec::<f64>::new(),
};

if args.is_ome {
if is_singlet {
raw = unravel(
ekore::operator_matrix_elements::unpolarized::spacelike::A_singlet(
args.order_qcd,
&mut c,
args.nf,
args.L,
),
args.order_qcd,
);
} else {
raw = unravel(
ekore::operator_matrix_elements::unpolarized::spacelike::A_non_singlet(
args.order_qcd,
&mut c,
args.nf,
args.L,
),
args.order_qcd,
);
}
} else if is_singlet {
raw = unravel(
ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd(
args.order_qcd,
&mut c,
args.nf,
),
args.order_qcd,
&mut c,
args.nf,
);
for gamma_s in res.iter().take(args.order_qcd) {
for col in gamma_s.iter().take(2) {
for el in col.iter().take(2) {
re.push(el.re);
im.push(el.im);
}
}
}
} else {
// we can not do 1D
let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qcd(
args.order_qcd,
args.mode0,
&mut c,
args.nf,
);
for el in res.iter().take(args.order_qcd) {
re.push(el.re);
im.push(el.im);
raw.re.push(el.re);
raw.im.push(el.im);
}
}

// pass on
(args.py)(
re.as_ptr(),
im.as_ptr(),
c.n.re,
c.n.im,
raw.re.as_ptr(),
raw.im.as_ptr(),
c.n().re,
c.n().im,
jac.re,
jac.im,
args.order_qcd,
Expand All @@ -72,6 +118,7 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 {
args.ev_op_max_order_qcd,
args.sv_mode_num,
args.is_threshold,
args.Lsv,
)
}

Expand Down Expand Up @@ -101,6 +148,7 @@ type PyQuadKerQCDT = unsafe extern "C" fn(
u8,
u8,
bool,
f64,
) -> f64;

/// Additional integration parameters
Expand Down Expand Up @@ -128,6 +176,8 @@ pub struct QuadQCDargs {
pub ev_op_max_order_qcd: u8,
pub sv_mode_num: u8,
pub is_threshold: bool,
pub is_ome: bool,
pub Lsv: f64,
}

/// Empty placeholder function for python callback.
Expand Down Expand Up @@ -159,14 +209,15 @@ pub unsafe extern "C" fn my_py(
_ev_op_max_order_qcd: u8,
_sv_mode_num: u8,
_is_threshold: bool,
_lsv: f64,
) -> f64 {
0.
}

/// Return empty additional arguments.
///
/// This is required to make the arguments part of the API, otherwise it won't be added to the compiled
/// package (since it does not appear in the signature of `quad_ker_qcd`).
/// package (since it does not appear in the signature of `rust_quad_ker_qcd`).
///
/// # Safety
/// This is the connection from and back to Python, so we don't know what is on the other side.
Expand All @@ -193,5 +244,7 @@ pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs {
ev_op_max_order_qcd: 0,
sv_mode_num: 0,
is_threshold: false,
is_ome: false,
Lsv: 0.,
}
}
13 changes: 13 additions & 0 deletions crates/ekore/refs.bib
Original file line number Diff line number Diff line change
Expand Up @@ -94,3 +94,16 @@ @article{Buza1996wv
pages = "301--320",
year = "1998"
}
@article{Bierenbaum2009zt,
author = "Bierenbaum, Isabella and Blumlein, Johannes and Klein, Sebastian",
title = "{The Gluonic Operator Matrix Elements at O(alpha(s)**2) for DIS Heavy Flavor Production}",
eprint = "0901.0669",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "DESY-08-187, SFB-CPP-08-107, IFIC-08-68",
doi = "10.1016/j.physletb.2009.01.057",
journal = "Phys. Lett. B",
volume = "672",
pages = "401--406",
year = "2009"
}
2 changes: 1 addition & 1 deletion crates/ekore/src/anomalous_dimensions.rs
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
//! The anomalous dimensions for DGLAP evolution.
//! The anomalous dimensions for |DGLAP| evolution.

pub mod unpolarized;
53 changes: 27 additions & 26 deletions crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
//! LO QCD
//! |LO| |QCD|.

use num::complex::Complex;

Expand All @@ -9,7 +9,7 @@ use crate::harmonics::cache::{Cache, K};
///
/// Implements Eq. (3.4) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
pub fn gamma_ns(c: &mut Cache, _nf: u8) -> Complex<f64> {
let N = c.n;
let N = c.n();
let S1 = c.get(K::S1);
let gamma = -(3.0 - 4.0 * S1 + 2.0 / N / (N + 1.0));
CF * gamma
Expand All @@ -19,7 +19,7 @@ pub fn gamma_ns(c: &mut Cache, _nf: u8) -> Complex<f64> {
///
/// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw).
pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
let N = c.n;
let N = c.n();
let gamma = -(N.powu(2) + N + 2.0) / (N * (N + 1.0) * (N + 2.0));
2.0 * TR * 2.0 * (nf as f64) * gamma
}
Expand All @@ -28,7 +28,7 @@ pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
///
/// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw).
pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex<f64> {
let N = c.n;
let N = c.n();
let gamma = -(N.powu(2) + N + 2.0) / (N * (N + 1.0) * (N - 1.0));
2.0 * CF * gamma
}
Expand All @@ -37,7 +37,7 @@ pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex<f64> {
///
/// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw).
pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
let N = c.n;
let N = c.n();
let S1 = c.get(K::S1);
let gamma = S1 - 1.0 / N / (N - 1.0) - 1.0 / (N + 1.0) / (N + 2.0);
CA * (4.0 * gamma - 11.0 / 3.0) + 4.0 / 3.0 * TR * (nf as f64)
Expand All @@ -54,63 +54,64 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {

#[cfg(test)]
mod tests {
use crate::cmplx;
use crate::{anomalous_dimensions::unpolarized::spacelike::as1::*, harmonics::cache::Cache};
use float_cmp::assert_approx_eq;
use super::*;
use crate::harmonics::cache::Cache;
use crate::{assert_approx_eq_cmplx, cmplx};
use num::complex::Complex;
use num::Zero;
const NF: u8 = 5;

#[test]
fn number_conservation() {
const N: Complex<f64> = cmplx![1., 0.];
const N: Complex<f64> = cmplx!(1., 0.);
let mut c = Cache::new(N);
let me = gamma_ns(&mut c, NF);
assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12);
assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12);
assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12);
}

#[test]
fn quark_momentum_conservation() {
const N: Complex<f64> = cmplx![2., 0.];
const N: Complex<f64> = cmplx!(2., 0.);
let mut c = Cache::new(N);
let me = gamma_ns(&mut c, NF) + gamma_gq(&mut c, NF);
assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12);
assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12);
assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12);
}

#[test]
fn gluon_momentum_conservation() {
const N: Complex<f64> = cmplx![2., 0.];
const N: Complex<f64> = cmplx!(2., 0.);
let mut c = Cache::new(N);
let me = gamma_qg(&mut c, NF) + gamma_gg(&mut c, NF);
assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12);
assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12);
assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12);
}

#[test]
fn gamma_qg_() {
const N: Complex<f64> = cmplx![1., 0.];
const N: Complex<f64> = cmplx!(1., 0.);
let mut c = Cache::new(N);
let me = gamma_qg(&mut c, NF);
assert_approx_eq!(f64, me.re, -20. / 3.0, ulps = 32);
assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12);
assert_approx_eq_cmplx!(f64, me, cmplx!(-20. / 3., 0.), ulps = 32, epsilon = 1e-12);
}

#[test]
fn gamma_gq_() {
const N: Complex<f64> = cmplx![0., 1.];
const N: Complex<f64> = cmplx!(0., 1.);
let mut c = Cache::new(N);
let me = gamma_gq(&mut c, NF);
assert_approx_eq!(f64, me.re, 4. / 3.0, ulps = 32);
assert_approx_eq!(f64, me.im, -4. / 3.0, ulps = 32);
assert_approx_eq_cmplx!(f64, me, cmplx!(4. / 3.0, -4. / 3.0), ulps = 32);
}

#[test]
fn gamma_gg_() {
const N: Complex<f64> = cmplx![0., 1.];
const N: Complex<f64> = cmplx!(0., 1.);
let mut c = Cache::new(N);
let me = gamma_gg(&mut c, NF);
assert_approx_eq!(f64, me.re, 5.195725159621, ulps = 32, epsilon = 1e-11);
assert_approx_eq!(f64, me.im, 10.52008856962, ulps = 32, epsilon = 1e-11);
assert_approx_eq_cmplx!(
f64,
me,
cmplx!(5.195725159621, 10.52008856962),
ulps = 32,
epsilon = 1e-11
);
}
}
Loading
Loading