+ 
+
+
+
+Please cite our DOI when using our code: 
## Contributing
Your feedback is welcome! If you want to report a (possible) bug or want to ask for a new feature, please raise an issue: 
\ No newline at end of file
diff --git a/doc/source/conf.py b/doc/source/conf.py
index f8b2c199f..a5eced7dc 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -23,9 +23,9 @@
author = 'N3PDF team'
# The short X.Y version
-version = '0.3'
+version = '0.4'
# The full version, including alpha/beta/rc tags
-release = 'v0.3.0'
+release = 'v0.4.0'
# -- General configuration ---------------------------------------------------
diff --git a/setup.py b/setup.py
index 23fe0b569..55f5e1e85 100644
--- a/setup.py
+++ b/setup.py
@@ -7,7 +7,7 @@
long_description = fh.read()
setup(name='eko',
- version='0.3.0',
+ version='0.4.0',
description='Evolution Kernel Operator',
long_description=long_description,
long_description_content_type="text/markdown",
diff --git a/src/eko/anomalous_dimensions/nlo.py b/src/eko/anomalous_dimensions/nlo.py
deleted file mode 100644
index 56c4e2fab..000000000
--- a/src/eko/anomalous_dimensions/nlo.py
+++ /dev/null
@@ -1,64 +0,0 @@
-# -*- coding: utf-8 -*-
-r"""
- This file contains the next-to-leading-order Altarelli-Parisi splitting kernels.
-"""
-
-#import numpy as np
-import numba as nb
-
-from eko import t_float, t_complex
-#from eko.ekomath import harmonic_S1 as S1
-
-
-@nb.njit
-def gamma_nsp_1(N: t_complex, nf: int, CA: t_float, CF: t_float):
- """
- Computes the next-to-leading-order non-singlet singlet-like anomalous dimension.
-
- Implements Eq. (3.5) of :cite:`Moch:2004pa`.
-
- Parameters
- ----------
- N : t_complex
- Mellin moment
- nf : int
- Number of active flavours
- CA : t_float
- Casimir constant of adjoint representation
- CF : t_float
- Casimir constant of fundamental representation
-
- Returns
- -------
- gamma_nsp_1 : t_complex
- Leading-order non-singlet anomalous dimension :math:`\\gamma_{ns}^{(1)+}(N)`
- """
- # TODO
- raise NotImplementedError("TODO")
-
-
-@nb.njit
-def gamma_nsm_1(N: t_complex, nf: int, CA: t_float, CF: t_float):
- """
- Computes the next-to-leading-order non-singlet valence-like anomalous dimension.
-
- Implements Eq. (3.6) of :cite:`Moch:2004pa`.
-
- Parameters
- ----------
- N : t_complex
- Mellin moment
- nf : int
- Number of active flavours
- CA : t_float
- Casimir constant of adjoint representation
- CF : t_float
- Casimir constant of fundamental representation
-
- Returns
- -------
- gamma_nsp_1 : t_complex
- Leading-order non-singlet anomalous dimension :math:`\\gamma_{ns}^{(1)-}(N)`
- """
- # TODO
- raise NotImplementedError("TODO")
diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py
index 309843870..c3c22ddca 100644
--- a/src/eko/strong_coupling.py
+++ b/src/eko/strong_coupling.py
@@ -7,6 +7,7 @@
"""
import numpy as np
+import scipy
import numba as nb
from eko import t_float
@@ -145,7 +146,7 @@ class StrongCoupling:
An instance of :class:`~eko.thresholds.ThresholdsConfig`
order: int
Evaluated order of the beta function: ``0`` = LO, ...
- method : ["analytic"]
+ method : ["expanded", "exact"]
Applied method to solve the beta function
Examples
@@ -161,7 +162,7 @@ class StrongCoupling:
"""
def __init__(
- self, consts, alpha_s_ref, scale_ref, thresh, order=0, method="analytic",
+ self, consts, alpha_s_ref, scale_ref, thresh, order=0, method="expanded",
):
# Sanity checks
if not isinstance(consts, constants.Constants):
@@ -175,7 +176,7 @@ def __init__(
if order not in [0, 1, 2]:
raise NotImplementedError("a_s beyond NNLO is not implemented")
self._order = order
- if method not in ["analytic"]:
+ if method not in ["expanded", "exact"]:
raise ValueError(f"Unknown method {method}")
self._method = method
self._constants = consts
@@ -226,9 +227,9 @@ def from_dict(cls, setup, constants, thresholds):
return cls(constants, alpha_ref, q2_alpha, thresholds, order)
# Hidden computation functions
- def _compute_analytic(self, as_ref, nf, scale_from, scale_to):
+ def _compute_expanded(self, as_ref, nf, scale_from, scale_to):
"""
- Compute via analytic expression.
+ Compute via expanded expression.
Parameters
----------
@@ -276,6 +277,51 @@ def _compute_analytic(self, as_ref, nf, scale_from, scale_to):
return res
+ def _compute_exact(self, as_ref, nf, scale_from, scale_to):
+ """
+ Compute via RGE.
+
+ Parameters
+ ----------
+ as_ref: t_float
+ reference alpha_s
+ nf: int
+ value of nf for computing alpha_s
+ scale_from: t_float
+ reference scale
+ scale_to : t_float
+ target scale
+
+ Returns
+ -------
+ a_s : t_float
+ strong coupling at target scale :math:`a_s(Q^2)`
+ """
+ # u = beta0 * ln(scale_to/scale_from)
+ beta0 = beta_0(nf, self._constants.CA, self._constants.CF, self._constants.TF)
+ u = beta0 * np.log(scale_to / scale_from)
+ b_vec = [1]
+ # NLO
+ if self._order >= 1:
+ beta1 = beta_1(
+ nf, self._constants.CA, self._constants.CF, self._constants.TF
+ )
+ b1 = beta1 / beta0
+ b_vec.append(b1)
+ # NNLO
+ if self._order >= 2:
+ beta2 = beta_2(
+ nf, self._constants.CA, self._constants.CF, self._constants.TF
+ )
+ b2 = beta2 / beta0
+ b_vec.append(b2)
+
+ def rge(_t, a, b_vec):
+ return -(a ** 2) * np.sum([a ** k * b for k, b in enumerate(b_vec)])
+
+ res = scipy.integrate.solve_ivp(rge, (0, u), (as_ref,), args=[b_vec])
+ return res.y[0][-1]
+
def _compute(self, as_ref, nf, scale_from, scale_to):
"""
Wrapper in order to pass the computation to the corresponding
@@ -298,8 +344,12 @@ def _compute(self, as_ref, nf, scale_from, scale_to):
strong coupling at target scale :math:`a_s(Q^2)`
"""
# TODO set up a cache system here
- # at the moment everything is analytic - and type has been checked in the constructor
- return self._compute_analytic(as_ref, nf, scale_from, scale_to)
+ # at the moment everything is expanded - and type has been checked in the constructor
+ if self._method == "exact":
+ as_new = self._compute_exact(as_ref, nf, scale_from, scale_to)
+ else:
+ as_new = self._compute_expanded(as_ref, nf, scale_from, scale_to)
+ return as_new
def a_s(self, scale_to, fact_scale=None):
"""
@@ -332,7 +382,7 @@ def a_s(self, scale_to, fact_scale=None):
if k < len(area_path) - 1:
next_nf_is_down = area_path[k + 1].nf < area.nf
# q2_to is the threshold value
- L = np.log(scale_to / fact_scale) # TODO why fact_scale instead of m2
+ L = np.log(scale_to / fact_scale)
if next_nf_is_down:
c1 = -4.0 / 3.0 * self._constants.TF * L
# TODO recover color constants
diff --git a/tests/test_ekomath.py b/tests/test_ekomath.py
index d5766c147..9bb73545a 100644
--- a/tests/test_ekomath.py
+++ b/tests/test_ekomath.py
@@ -1,5 +1,6 @@
# -*- coding: utf-8 -*-
from numpy.testing import assert_almost_equal
+import pytest
from eko import t_complex
from eko import ekomath
@@ -73,6 +74,12 @@ def test_cern_polygamma():
me = ekomath.cern_polygamma(z, k)
ref = fortran_ref[nk][nz]
assert_almost_equal(me, ref)
+ # errors
+ with pytest.raises(NotImplementedError):
+ _ = ekomath.cern_polygamma(1, 5)
+ with pytest.raises(ValueError):
+ _ = ekomath.cern_polygamma(0, 0)
+
def test_harmonic_S1():
diff --git a/tests/test_strong_coupling.py b/tests/test_strong_coupling.py
index f09bbfc00..6fb802236 100644
--- a/tests/test_strong_coupling.py
+++ b/tests/test_strong_coupling.py
@@ -10,6 +10,7 @@
from eko import thresholds
from eko.constants import Constants
+# try to load LHAPDF - if not available, we'll use the cached values
try:
import lhapdf
@@ -17,6 +18,7 @@
except ImportError:
use_LHAPDF = False
+# try to load APFEL - if not available, we'll use the cached values
try:
import apfel
@@ -131,6 +133,30 @@ def test_ref(self):
sc.a_s(scale_ref), alphas_ref / 4.0 / np.pi
)
+ def test_exact_LO(self):
+ # prepare
+ thresh_setups = [
+ {"FNS": "FFNS", "NfFF": 3},
+ {"FNS": "FFNS", "NfFF": 4},
+ {"FNS": "ZM-VFNS", "mc": 2, "mb": 4, "mt": 175},
+ ]
+ alphas_ref = 0.118
+ scale_ref = 91.0 ** 2
+ for thresh_setup in thresh_setups:
+ thresh_setup["Q0"] = 1
+ thresholds_conf = thresholds.ThresholdsConfig.from_dict(thresh_setup)
+ # in LO expanded = exact
+ sc_expanded = StrongCoupling(
+ constants, alphas_ref, scale_ref, thresholds_conf, 0, "expanded"
+ )
+ sc_exact = StrongCoupling(
+ constants, alphas_ref, scale_ref, thresholds_conf, 0, "exact"
+ )
+ for q2 in [1, 1e1, 1e2, 1e3, 1e4]:
+ np.testing.assert_allclose(
+ sc_expanded.a_s(q2), sc_exact.a_s(q2), rtol=5e-4
+ )
+
class BenchmarkStrongCoupling:
def test_a_s(self):
@@ -357,6 +383,135 @@ def benchmark_lhapdf_ffns_lo(self):
# check myself to LHAPDF
np.testing.assert_allclose(lhapdf_vals, np.array(my_vals), rtol=5e-4)
+ def benchmark_apfel_exact(self):
+ """test exact towards APFEL"""
+ Q2s = [1e1, 1e2, 1e3, 1e4]
+ alphas_ref = 0.118
+ scale_ref = 90 ** 2
+ # collect my values
+ threshold_holder = thresholds.ThresholdsConfig(scale_ref, "FFNS", nf=3)
+ # LHAPDF cache
+ apfel_vals_dict = {
+ 0: np.array(
+ [
+ 0.021635019899707245,
+ 0.014937719308417242,
+ 0.01140668497649318,
+ 0.009225844999427163,
+ ]
+ ),
+ 1: np.array(
+ [
+ 0.025027723843261098,
+ 0.015730685887616093,
+ 0.01159096381106341,
+ 0.009215179564010682,
+ ]
+ ),
+ 2: np.array(
+ [
+ 0.025717091835015565,
+ 0.01583723253352162,
+ 0.011610857909393214,
+ 0.009214183434685514,
+ ]
+ ),
+ }
+ for order in range(2 + 1):
+ sc = StrongCoupling(
+ constants,
+ alphas_ref,
+ scale_ref,
+ threshold_holder,
+ order=order,
+ method="exact",
+ )
+ my_vals = []
+ for Q2 in Q2s:
+ my_vals.append(sc.a_s(Q2))
+ # get APFEL numbers - if available else use cache
+ apfel_vals = apfel_vals_dict[order]
+ if use_APFEL:
+ # run apfel
+ apfel.CleanUp()
+ apfel.SetTheory("QCD")
+ apfel.SetPerturbativeOrder(order)
+ apfel.SetAlphaEvolution("exact")
+ apfel.SetAlphaQCDRef(alphas_ref, np.sqrt(scale_ref))
+ apfel.SetFFNS(3)
+ apfel.SetRenFacRatio(1)
+ apfel.InitializeAPFEL()
+ # collect a_s
+ apfel_vals_cur = []
+ for Q2 in Q2s:
+ apfel_vals_cur.append(apfel.AlphaQCD(np.sqrt(Q2)) / (4.0 * np.pi))
+ # print(apfel_vals_cur)
+ np.testing.assert_allclose(apfel_vals, np.array(apfel_vals_cur))
+ # check myself to APFEL
+ np.testing.assert_allclose(apfel_vals, np.array(my_vals), rtol=2e-4)
+
+ def benchmark_lhapdf_exact(self):
+ """test exact towards LHAPDF"""
+ Q2s = [1e1, 1e2, 1e3, 1e4]
+ alphas_ref = 0.118
+ scale_ref = 90 ** 2
+ # collect my values
+ threshold_holder = thresholds.ThresholdsConfig(scale_ref, "FFNS", nf=3)
+ # LHAPDF cache
+ lhapdf_vals_dict = {
+ 0: np.array(
+ [
+ 0.021634715590772086,
+ 0.014937700253543466,
+ 0.011406683936848657,
+ 0.009225845071914008,
+ ]
+ ),
+ 1: np.array(
+ [
+ 0.025027589743439077,
+ 0.015730666136188308,
+ 0.011590962671090168,
+ 0.009215179641099218,
+ ]
+ ),
+ 2: np.array(
+ [
+ 0.02571700975829869,
+ 0.015837212322787945,
+ 0.011610856755395813,
+ 0.009214183512196827,
+ ]
+ ),
+ }
+ for order in range(2 + 1):
+ sc = StrongCoupling(
+ constants,
+ alphas_ref,
+ scale_ref,
+ threshold_holder,
+ order=order,
+ method="exact",
+ )
+ my_vals = []
+ for Q2 in Q2s:
+ my_vals.append(sc.a_s(Q2))
+ lhapdf_vals = lhapdf_vals_dict[order]
+ if use_LHAPDF:
+ as_lhapdf = lhapdf.mkBareAlphaS("ODE")
+ as_lhapdf.setOrderQCD(1 + order)
+ as_lhapdf.setFlavorScheme("FIXED", 3)
+ as_lhapdf.setAlphaSMZ(alphas_ref)
+ as_lhapdf.setMZ(np.sqrt(scale_ref))
+ # collect a_s
+ lhapdf_vals_cur = []
+ for Q2 in Q2s:
+ lhapdf_vals_cur.append(as_lhapdf.alphasQ2(Q2) / (4.0 * np.pi))
+ # print(lhapdf_vals_cur)
+ np.testing.assert_allclose(lhapdf_vals, np.array(lhapdf_vals_cur))
+ # check
+ np.testing.assert_allclose(lhapdf_vals, np.array(my_vals), rtol=2e-4)
+
def benchmark_lhapdf_zmvfns_lo(self):
"""test ZM-VFNS LO towards LHAPDF"""
Q2s = [1, 1e1, 1e2, 1e3, 1e4]
@@ -369,13 +524,13 @@ def benchmark_lhapdf_zmvfns_lo(self):
# compute all Lambdas
# Lambda2_5 = self._get_Lambda2_LO(alphas_ref / (4.0 * np.pi), scale_ref, 5)
# as_FFNS_LO_5 = StrongCoupling(
- # constants, alphas_ref, scale_ref, 0, "FFNS", nf=5, method="analytic"
+ # constants, alphas_ref, scale_ref, 0, "FFNS", nf=5, method="expanded"
# )
# Lambda2_6 = self._get_Lambda2_LO(as_FFNS_LO_5.a_s(m2t), m2t, 6)
# as_b = as_FFNS_LO_5.a_s(m2b)
# Lambda2_4 = self._get_Lambda2_LO(as_b, m2b, 4)
# as_FFNS_LO_4 = StrongCoupling(
- # constants, as_b * 4.0 * np.pi, m2b, 0, "FFNS", nf=4, method="analytic"
+ # constants, as_b * 4.0 * np.pi, m2b, 0, "FFNS", nf=4, method="expanded"
# )
# Lambda2_3 = self._get_Lambda2_LO(as_FFNS_LO_4.a_s(m2c), m2c, 3)