Merge branch 'master' into rust

.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

@@ -5,7 +5,7 @@ ci:
   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
@@ -14,16 +14,16 @@ repos:
       - id: check-merge-conflict
       - id: debug-statements
   - repo: https://github.com/hadialqattan/pycln
-    rev: v2.2.2
+    rev: v2.4.0
     hooks:
       - id: pycln
         args: [--config=pyproject.toml]
   - repo: https://github.com/psf/black
-    rev: 23.7.0
+    rev: 23.11.0
     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
@@ -32,7 +32,7 @@ repos:
       - id: isort
         args: ["--profile", "black"]
   - repo: https://github.com/asottile/pyupgrade
-    rev: v3.10.1
+    rev: v3.15.0
     hooks:
       - id: pyupgrade
   - repo: https://github.com/pycqa/pydocstyle
@@ -62,6 +62,6 @@ repos:
       files: ^crates/ekore/.*\.rs$
       args: []
   - repo: https://github.com/pre-commit/pre-commit
-    rev: v3.3.3
+    rev: v3.5.0
     hooks:
       - id: validate_manifest

benchmarks/eko/benchmark_inverse_matching.py

@@ -0,0 +1,107 @@
+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],
+)
+
+# 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]
+        )

doc/source/refs.bib

@@ -999,3 +999,36 @@ @article{Falcioni:2023vqq
     month = "7",
     year = "2023"
 }
+
+@article{Gehrmann:2023cqm,
+    author = "Gehrmann, Thomas and von Manteuffel, Andreas and Sotnikov, Vasily and Yang, Tong-Zhi",
+    title = "{Complete $N_f^2$ contributions to four-loop pure-singlet splitting functions}",
+    eprint = "2308.07958",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "MSUHEP-23-024, ZU-TH 43/23",
+    month = "8",
+    year = "2023"
+}
+
+@article{Falcioni:2023tzp,
+    author = "Falcioni, G. and Herzog, F. and Moch, S. and Vermaseren, J. and Vogt, A.",
+    title = "{The double fermionic contribution to the four-loop quark-to-gluon splitting function}",
+    eprint = "2310.01245",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "ZU-TH 62/23, DESY 23-146, Nikhef 2023-015, LTH 1353",
+    month = "10",
+    year = "2023"
+}
+
+@article{Moch:2023tdj,
+    author = "Moch, S. and Ruijl, B. and Ueda, T. and Vermaseren, J. and Vogt, A.",
+    title = "{Additional moments and x-space approximations of four-loop splitting functions in QCD}",
+    eprint = "2310.05744",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "DESY-23-150, Nikhef 23-016, LTH 1354",
+    month = "10",
+    year = "2023"
+}

doc/source/theory/Mellin.rst

@@ -56,7 +56,7 @@ Plus Distributions
 ^^^^^^^^^^^^^^^^^^
 
 .. math ::
-    \mathcal M[1/(1-x)_+](N) = S_1(N)
+    \mathcal M[1/(1-x)_+](N) = - S_1(N - 1)
 
 with the harmonic sum :math:`S_1` (see :ref:`theory/mellin:harmonic sums`).
 

doc/source/theory/N3LO_ad.rst

@@ -90,16 +90,16 @@ In |EKO| they are implemented as follows:
         -   The large-N limit :cite:`Moch:2017uml`, which reads (Eq. 2.17):
 
             .. math ::
-                \gamma_{ns} \approx A^{(f)}_4 S_1(N) - B_4 + C_4 \frac{S_1(N)}{N} - (D_4 + \frac{1}{2} A^{(f)}_4) \frac{1}{N} + \mathcal{O}(\frac{\ln^k(N)}{N^2})
+                \gamma_{ns} \approx A^{(f)}_4 S_1(N) - B^{(f)}_4 + C^{(f)}_4 \frac{S_1(N)}{N} - D^{(f)}_4 \frac{1}{N}
 
             This limit is common for all :math:`\gamma_{ns,+}^{(3)},\gamma_{ns,-}^{(3)},\gamma_{ns,v}^{(3)}`.
             The coefficient :math:`A^{(f)}_4`, being related to the twist-2 spin-N operators,
             can be obtained from the |QCD| cusp calculation
-            :cite:`Henn:2019swt`, while the :math:`B_4` is fixed by the integral of the 4-loop splitting function
+            :cite:`Henn:2019swt`, while the :math:`B^{(f)}_4` is fixed by the integral of the 4-loop splitting function
             and has been firstly computed in :cite:`Moch:2017uml` in the large :math:`n_c` limit.
             More recently :cite:`Duhr:2022cob`, it has been determined  in the full color expansion
             by computing various |N3LO| cross sections in the soft limit.
-            :math:`C_4,D_4` instead can be computed directly from lower order splitting functions.
+            :math:`C^{(f)}_4,D^{(f)}_4` instead can be computed directly from lower order splitting functions.
             From large-x resummation :cite:`Davies:2016jie`, it is possible to infer further constrains
             on sub-leading terms :math:`\frac{\ln^k(N)}{N^2}`, since the non-singlet splitting
             functions contain only terms :math:`(1-x)^a\ln^k(1-x)` with :math:`a \ge 1`.
@@ -132,7 +132,7 @@ In |EKO| they are implemented as follows:
                 *   - :math:`\frac{x}{2}\ln^2(x)`
                     - :math:`\frac{1}{(N+1)^3}`
                 *   - :math:`x^{2}, x^{3}`
-                    - :math:`\frac{1}{(N-2)},\frac{1}{(N-3)}`
+                    - :math:`\frac{1}{(N+2)},\frac{1}{(N+3)}`
 
             The first five functions model the sub-leading differences in the :math:`N\to \infty` limit,
             while the last three help the convergence in the small-N region. Finally, we add a polynomial part
@@ -193,8 +193,13 @@ the following terms:
             - |T|
             - |T|
 
-Only the parts proportional to :math:`n_f^3` are known analytically
+The parts proportional to :math:`n_f^3` are known analytically
 :cite:`Davies:2016jie` and have been included so far.
+For :math:`\gamma_{qq,ps}` and :math:`\gamma_{gq}` also the component
+proportional to :math:`n_f^2` has been computed in :cite:`Gehrmann:2023cqm`
+and :cite:`Falcioni:2023tzp` respectively and it's used in our code
+through an approximations obtained with 30 moments.
+
 The other parts are approximated using some known limits:
 
     *   The small-x limit, given in the large :math:`N_c` approximation by
@@ -229,11 +234,11 @@ The other parts are approximated using some known limits:
         It is known that :cite:`Albino:2000cp,Moch:2021qrk` the diagonal terms diverge in N-space as:
 
             .. math ::
-                \gamma_{kk} \approx A^{(r)}_4 S_1(N) + B^{(r)}_4 + C^{(r)}_4 \frac{S_1(N)}{N} + \mathcal{O}(\frac{1}{N})
+                \gamma_{kk} \approx A^{(r)}_4 S_1(N) + B^{(r)}_4 + C^{(r)}_4 \frac{S_1(N)}{N} - D^{(r)}_4 \frac{1}{N}
 
         Where again the coefficient :math:`A^{(r)}_4` is the |QCD| cusp anomalous dimension for the adjoint or fundamental representation,
         the coefficient :math:`B^{(r)}_4` has been extracted from soft anomalous dimensions :cite:`Duhr:2022cob`.
-        and :math:`C^{(r)}_4`can be estimate from lower orders :cite:`Dokshitzer:2005bf`.
+        and :math:`C^{(r)}_4,D^{(r)}_4`can be estimate from lower orders :cite:`Dokshitzer:2005bf`.
         However, :math:`\gamma_{qq,ps}^{(3)}` do not constrain any divergence at large-x or constant term so its expansion starts as
         :math:`\mathcal{O}(\frac{1}{N^2})`.
         The off-diagonal do not contain any +-distributions or delta distributions but can include divergent logarithms
@@ -252,14 +257,14 @@ The other parts are approximated using some known limits:
                 \gamma_{qq,ps} \approx (1-x)[c_{4} \ln^4(1-x) + c_{3} \ln^3(1-x)] + \mathcal{O}((1-x)\ln^2(1-x))
 
 
-    *   The 4 lowest even N moments provided in :cite:`Moch:2021qrk`, where we can use momentum conservation
-        to fix:
+    *   The 5 lowest even N moments provided in :cite:`Moch:2021qrk,Moch:2023tdj`,
+        where momentum conservation fixes:
 
             .. math ::
                 & \gamma_{qg}(2) + \gamma_{gg}(2) = 0 \\
                 & \gamma_{qq}(2) + \gamma_{gq}(2) = 0 \\
 
-        For :math:`\gamma_{qq,ps}, \gamma_{qg}` other 6 additional moments are available :cite:`Falcioni:2023luc,Falcioni:2023vqq`.
+        For :math:`\gamma_{qq,ps}, \gamma_{qg}` other 5 additional moments are available :cite:`Falcioni:2023luc,Falcioni:2023vqq`.
         making the parametrization of this splitting function much more accurate.
 
 The difference between the known moments and the known limits is parametrized
@@ -276,9 +281,9 @@ we need to account for a possible source of uncertainties arising during the app
 This uncertainty is neglected in the non-singlet case.
 
 The procedure is performed in two steps for each different anomalous dimension separately.
-First, we solve the system associated to the 4 known moments,
+First, we solve the system associated to the 5 (10) known moments,
 minus the known limits, using different functional bases.
-Any possible candidate contains 4 elements and is obtained with the following prescription:
+Any possible candidate contains 5 elements and is obtained with the following prescription:
 
     1. one function is leading small-N unknown contribution, which correspond to the highest power unknown for the pole at :math:`N=1`,
 
@@ -312,32 +317,29 @@ final reduced sets of candidates.
         :align: center
 
         *   - :math:`f_1(N)`
-            - :math:`\frac{S_2(N-2)}{N}`
+            - :math:`\frac{1}{(N-1)^2}`
         *   - :math:`f_2(N)`
-            - :math:`\frac{1}{N}`
+            - :math:`\mathcal{M}[(1-x)\ln^3(1-x)]`
         *   - :math:`f_3(N)`
-            - :math:`\frac{1}{N-1},\ \frac{S_1(N)}{N^2}`
+            - :math:`\frac{1}{N-1},`
         *   - :math:`f_4(N)`
-            - :math:`\frac{1}{N-1},\ \frac{1}{N^4},\ \frac{1}{N^3},\ \frac{1}{N^2},\ \frac{1}{(N+1)^3},\ \frac{1}{(N+1)^2},\ \frac{1}{N+1},\ \frac{1}{N+2},\ \mathcal{M}[(1-x)\ln(1-x)],\ \frac{S_1(N)}{N^2}, \ \mathcal{M}[(1-x)^2\ln(1-x)],`
+            - :math:`\frac{1}{N^4},\ \frac{1}{N^3},\ \frac{1}{N^2},\ \frac{1}{(N+1)},\ \frac{1}{(N+2)},\ \mathcal{M}[(1-x)\ln^2(1-x)],\ \mathcal{M}[(1-x)\ln(1-x)]`
 
     .. list-table::  :math:`\gamma_{gq}^{(3)}` parametrization basis
         :align: center
 
         *   - :math:`f_1(N)`
-            - :math:`\frac{S_2(N-2)}{N}`
+            - :math:`\frac{1}{(N-1)^2}`
         *   - :math:`f_2(N)`
-            - :math:`\frac{S_1^3(N)}{N}`
+            - :math:`\mathcal{M}[\ln^3(1-x)]`
         *   - :math:`f_3(N)`
-            - :math:`\frac{1}{N-1},\ \frac{1}{N^4}`
+            - :math:`\frac{1}{N-1}`
         *   - :math:`f_4(N)`
-            - :math:`\frac{1}{N-1},\ \frac{1}{N^4},\ \frac{1}{N^3},\ \frac{1}{N^2},\ \frac{1}{N},\ \frac{1}{(N+1)^3},\ \frac{1}{(N+1)^2},\ \frac{1}{N+1},\ \frac{1}{N+2},\ \frac{S_1(N-2)}{N},\ \mathcal{M}[\ln^3(1-x)],\ \mathcal{M}[\ln^2(1-x)], \frac{S_1(N)}{N},\ \frac{S_1^2(N)}{N}`
-
-    Note that this table refers only to the :math:`n_f^0` part where we assume no violation of the scaling with :math:`\gamma_{gg}`
-    also for the |NLL| term, to help the convergence. We expect that any possible deviation can be parametrized as a shift in the |NNLL| terms
-    and in the |NLL| :math:`n_f^1` which are free to vary independently.
-    Furthermore for the part :math:`\propto n_f^2` we adopt a slightly different
-    basis to account fot the fact that the leading
-    contribution for the pole at :math:`N=1` is :math:`\frac{1}{(N-1)^2}`.
+            - :math:`\frac{1}{N^4},\ \frac{1}{N^3},\ \frac{1}{N^2},\ \frac{1}{(N+1)},\ \frac{1}{(N+2)},\ \mathcal{M}[\ln^2(1-x)],\ \mathcal{M}[\ln(1-x)]`
+
+    Following :cite:`Moch:2023tdj` we have assumed no violation of the scaling with :math:`\gamma_{gg}`
+    also for the |NLL| small-x term, to help the convergence. We expect that any possible deviation can be parametrized as a shift in the |NNLL| terms
+    which are free to vary independently.
 
 Slightly different choices are performed for :math:`\gamma_{gq}^{(3)}` and :math:`\gamma_{qq,ps}^{(3)}`
 where 10 moments are known. In this case we can select a larger number of functions in group 3

extras/lh_bench_23/.gitignore

@@ -0,0 +1,2 @@
+*.tar
+*.csv