Merge branch 'master' into load-opened

.envrc

@@ -0,0 +1,11 @@
+if ! has nix_direnv_version || ! nix_direnv_version 2.2.1; then
+  source_url "https://raw.githubusercontent.com/nix-community/nix-direnv/2.2.1/direnvrc" "sha256-zelF0vLbEl5uaqrfIzbgNzJWGmLzCmYAkInj/LNxvKs="
+fi
+
+nix_direnv_watch_file devenv.nix
+nix_direnv_watch_file devenv.lock
+nix_direnv_watch_file devenv.yaml
+if ! use flake . --impure
+then
+  echo "devenv could not be built. The devenv environment was not loaded. Make the necessary changes to devenv.nix and hit enter to try again." >&2
+fi

.gitignore

@@ -160,3 +160,4 @@ validation-report.json
 
 # Folders to ignore
 node_modules
+.devenv

.pre-commit-config.yaml

@@ -2,9 +2,10 @@
 # See https://pre-commit.com/hooks.html for more hooks
 ci:
   autofix_prs: false
+  skip: [fmt-eko, fmt-ekore]
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.4.0
+    rev: v4.5.0
     hooks:
       - id: trailing-whitespace
       - id: end-of-file-fixer
@@ -13,25 +14,25 @@ repos:
       - id: check-merge-conflict
       - id: debug-statements
   - repo: https://github.com/hadialqattan/pycln
-    rev: v2.1.6
+    rev: v2.4.0
     hooks:
       - id: pycln
         args: [--config=pyproject.toml]
-  - repo: https://github.com/psf/black
-    rev: 23.7.0
+  - repo: https://github.com/psf/black-pre-commit-mirror
+    rev: 23.12.1
     hooks:
       - id: black
   - repo: https://github.com/asottile/blacken-docs
-    rev: 1.15.0
+    rev: 1.16.0
     hooks:
       - id: blacken-docs
   - repo: https://github.com/pycqa/isort
-    rev: 5.12.0
+    rev: 5.13.2
     hooks:
       - id: isort
         args: ["--profile", "black"]
   - repo: https://github.com/asottile/pyupgrade
-    rev: v3.9.0
+    rev: v3.15.0
     hooks:
       - id: pyupgrade
   - repo: https://github.com/pycqa/pydocstyle
@@ -42,7 +43,25 @@ repos:
         args: ["--add-ignore=D107,D105"]
         additional_dependencies:
           - toml
+  - repo: local
+    hooks:
+    - id: fmt-eko
+      name: fmt-eko
+      description: Format eko files with cargo fmt.
+      entry: cargo fmt --manifest-path crates/eko/Cargo.toml --
+      language: system
+      files: ^crates/eko/.*\.rs$
+      args: []
+  - repo: local
+    hooks:
+    - id: fmt-ekore
+      name: fmt-ekore
+      description: Format ekore files with cargo fmt.
+      entry: cargo fmt --manifest-path crates/ekore/Cargo.toml --
+      language: system
+      files: ^crates/ekore/.*\.rs$
+      args: []
   - repo: https://github.com/pre-commit/pre-commit
-    rev: v3.3.3
+    rev: v3.6.0
     hooks:
       - id: validate_manifest

benchmarks/eko/benchmark_alphaem.py

@@ -0,0 +1,144 @@
+"""This module benchmarks alpha_em against alphaQED23 and validphys.
+
+alphaQED23 can be obtained from http://www-com.physik.hu-berlin.de/~fjeger/software.html .
+"""
+
+import numpy as np
+import pytest
+
+from eko.couplings import Couplings
+from eko.quantities.couplings import CouplingEvolutionMethod, CouplingsInfo
+from eko.quantities.heavy_quarks import QuarkMassScheme
+
+
+@pytest.mark.isolated
+class BenchmarkCouplings:
+    def test_alphaQED_high(self):
+        """testing aginst alphaQED23 for high Q values"""
+        alphaQED23 = np.array(
+            [
+                0.0075683,
+                0.0075867,
+                0.0076054,
+                0.0076245,
+                0.0076437,
+                0.0076631,
+                0.0076827,
+                0.0077024,
+                0.0077222,
+                0.0077421,
+                0.0077621,
+            ]
+        )
+        qvec = np.geomspace(10, 100, 11)
+        couplings = CouplingsInfo.from_dict(
+            dict(
+                alphas=0.118,
+                alphaem=7.7553e-03,
+                scale=91.2,
+                num_flavs_ref=5,
+                max_num_flavs=5,
+                em_running=True,
+            )
+        )
+        eko_alpha = Couplings(
+            couplings,
+            (3, 2),
+            method=CouplingEvolutionMethod.EXACT,
+            masses=[m**2 for m in [1.51, 4.92, 172.5]],
+            hqm_scheme=QuarkMassScheme.POLE,
+            thresholds_ratios=[1.0, 1.0, np.inf],
+        )
+        alpha_eko = np.array([eko_alpha.a_em(q**2) * 4 * np.pi for q in qvec])
+
+        np.testing.assert_allclose(alphaQED23, alpha_eko, rtol=1.8e-4)
+
+    def test_alphaQED_low(self):
+        """testing aginst alphaQED23 for low Q values: they are close but not identical"""
+        alphaQED23 = np.array(
+            [
+                0.0074192,
+                0.007431,
+                0.0074434,
+                0.0074565,
+                0.0074702,
+                0.0074847,
+                0.0075,
+                0.0075161,
+                0.0075329,
+                0.0075503,
+                0.0075683,
+            ]
+        )
+        qvec = np.geomspace(1, 10, 11)
+        couplings = CouplingsInfo.from_dict(
+            dict(
+                alphas=0.118,
+                alphaem=7.7553e-03,
+                scale=91.2,
+                num_flavs_ref=5,
+                max_num_flavs=5,
+                em_running=True,
+            )
+        )
+        eko_alpha = Couplings(
+            couplings,
+            (3, 2),
+            method=CouplingEvolutionMethod.EXACT,
+            masses=[m**2 for m in [1.51, 4.92, 172.5]],
+            hqm_scheme=QuarkMassScheme.POLE,
+            thresholds_ratios=[1.0, 1.0, np.inf],
+        )
+        alpha_eko = np.array([eko_alpha.a_em(q**2) * 4 * np.pi for q in qvec])
+
+        np.testing.assert_allclose(alphaQED23, alpha_eko, rtol=3.2e-3)
+
+    def test_validphys(self):
+        """testing aginst validphys"""
+        alpha_vp = np.array(
+            [
+                0.007408255774054356,
+                0.007425240094018394,
+                0.007449051094996458,
+                0.007476301027742958,
+                0.007503751810862984,
+                0.007532299008699658,
+                0.007561621958607614,
+                0.007591174885612722,
+                0.007620960508577136,
+                0.00765098158940323,
+                0.007681240933888789,
+                0.007711741392602655,
+                0.007742485861781425,
+                0.007773477284247778,
+                0.007804718650351058,
+                0.007836212998930739,
+                0.00786796341830342,
+                0.007899973047274033,
+                0.007932245076171957,
+                0.00796478274791273,
+            ]
+        )
+        qvec = np.geomspace(1, 1000, 20)
+        couplings = CouplingsInfo.from_dict(
+            dict(
+                alphas=0.118,
+                alphaem=7.7553e-03,
+                scale=91.2,
+                num_flavs_ref=5,
+                max_num_flavs=5,
+                em_running=True,
+            )
+        )
+        eko_alpha = Couplings(
+            couplings,
+            (3, 2),
+            method=CouplingEvolutionMethod.EXACT,
+            masses=[m**2 for m in [1.51, 4.92, 172.5]],
+            hqm_scheme=QuarkMassScheme.POLE,
+            thresholds_ratios=[1.0, 1.0, np.inf],
+        )
+        eko_alpha.decoupled_running = True
+        alpha_eko = np.array([eko_alpha.a_em(q**2) * 4 * np.pi for q in qvec])
+
+        np.testing.assert_allclose(alpha_vp, alpha_eko, rtol=5e-6)

benchmarks/eko/benchmark_inverse_matching.py

@@ -0,0 +1,108 @@
+import pathlib
+
+import numpy as np
+import pytest
+from banana import toy
+
+from eko import EKO
+from eko.io import runcards
+from eko.io.types import ReferenceRunning
+from eko.runner.managed import solve
+from ekobox import apply
+
+here = pathlib.Path(__file__).parent.absolute()
+MC = 1.51
+C_PID = 4
+
+# theory settings
+th_raw = dict(
+    order=[3, 0],
+    couplings=dict(
+        alphas=0.118,
+        alphaem=0.007496252,
+        scale=91.2,
+        num_flavs_ref=5,
+        max_num_flavs=6,
+    ),
+    heavy=dict(
+        num_flavs_init=4,
+        num_flavs_max_pdf=6,
+        intrinsic_flavors=[],
+        masses=[ReferenceRunning([mq, np.nan]) for mq in (MC, 4.92, 172.5)],
+        masses_scheme="POLE",
+        matching_ratios=[1.0, 1.0, np.inf],
+    ),
+    xif=1.0,
+    n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0),
+    matching_order=[2, 0],
+    use_fhmruvv=False,
+)
+
+# operator settings
+op_raw = dict(
+    mu0=1.65,
+    xgrid=[0.0001, 0.001, 0.01, 0.1, 1],
+    mugrid=[(MC, 3), (MC, 4)],
+    configs=dict(
+        evolution_method="truncated",
+        ev_op_max_order=[1, 0],
+        ev_op_iterations=1,
+        interpolation_polynomial_degree=4,
+        interpolation_is_log=True,
+        scvar_method="exponentiated",
+        inversion_method="exact",
+        n_integration_cores=0,
+        polarized=False,
+        time_like=False,
+    ),
+    debug=dict(
+        skip_singlet=False,
+        skip_non_singlet=False,
+    ),
+)
+
+
+@pytest.mark.isolated
+def benchmark_inverse_matching():
+    th_card = runcards.TheoryCard.from_dict(th_raw)
+    op_card = runcards.OperatorCard.from_dict(op_raw)
+
+    eko_path2 = here / "test2.tar"
+    eko_path2.unlink(missing_ok=True)
+    solve(th_card, op_card, eko_path2)
+
+    th_card.matching_order = [1, 0]
+    eko_path1 = here / "test1.tar"
+    eko_path1.unlink(missing_ok=True)
+    solve(th_card, op_card, eko_path1)
+
+    eko_output1 = EKO.read(eko_path1)
+    eko_output2 = EKO.read(eko_path2)
+    op1_nf3 = eko_output1[(MC**2, 3)]
+    op2_nf3 = eko_output2[(MC**2, 3)]
+    op1_nf4 = eko_output1[(MC**2, 4)]
+    op2_nf4 = eko_output2[(MC**2, 4)]
+
+    # test that nf=4 operators are the same
+    np.testing.assert_allclose(op1_nf4.operator, op2_nf4.operator)
+
+    with pytest.raises(AssertionError):
+        np.testing.assert_allclose(op2_nf3.operator, op2_nf4.operator)
+
+    with pytest.raises(AssertionError):
+        np.testing.assert_allclose(op1_nf3.operator, op1_nf4.operator)
+
+    with pytest.raises(AssertionError):
+        np.testing.assert_allclose(op1_nf3.operator, op2_nf3.operator)
+
+    pdf1 = apply.apply_pdf(eko_output1, toy.mkPDF("ToyLH", 0))
+    pdf2 = apply.apply_pdf(eko_output2, toy.mkPDF("ToyLH", 0))
+
+    # test that different PTO matching is applied correctly
+    np.testing.assert_allclose(
+        pdf1[(MC**2, 4)]["pdfs"][C_PID], pdf2[(MC**2, 4)]["pdfs"][C_PID]
+    )
+    with pytest.raises(AssertionError):
+        np.testing.assert_allclose(
+            pdf1[(MC**2, 3)]["pdfs"][C_PID], pdf2[(MC**2, 3)]["pdfs"][C_PID]
+        )

crates/.gitignore

@@ -0,0 +1 @@
+/target