N3LO matching

benchmarks/.gitignore

@@ -1,3 +1,7 @@
+# ignore airspeed velocity output
+html/
+results/
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]

benchmarks/NNPDF_bench.py

@@ -1,9 +1,14 @@
 # -*- coding: utf-8 -*-
 """
-Benchmark NNDPF31 pdf family
+Benchmark NNPDF pdf family
 """
+import numpy as np
+from banana import register
+
 from ekomark.benchmark.runner import Runner
 
+register(__file__)
+
 
 class BenchmarkNNPDF(Runner):
     """
@@ -26,9 +31,7 @@ def skip_pdfs(self, _theory):
         ]
 
 
-base_operator = {
-    "ev_op_iterations": 1,
-}
+base_operator = {"ev_op_iterations": 1, "backward_inversion": "exact"}
 
 base_theory = {
     "Qref": 91.2,
@@ -39,7 +42,6 @@ def skip_pdfs(self, _theory):
     "kbThr": 1.0,
     "ktThr": 1.0,
     "alphas": 0.118000,
-    "FNS": "ZM-VFNS",
     "ModEv": "TRN",
 }
 
@@ -75,6 +77,8 @@ def benchmark_nnlo(self, Q0=1.65, Q2grid=(100,)):
         theory_card = {
             **base_theory,
             "PTO": 2,
+            "IC": 1,
+            "IB": 1,
             "Q0": Q0,
         }
 
@@ -83,13 +87,13 @@ def benchmark_nnlo(self, Q0=1.65, Q2grid=(100,)):
 
 
 if __name__ == "__main__":
+    low2 = 4**2
+    high2 = 30**2
     # nn31 = BenchmarkNNPDF31()
-    # low2 = 1.65**2
-    # high2 = 20
     # # test forward
     # nn31.benchmark_nlo(Q0=np.sqrt(low2), Q2grid=[10])
     # # test backward
     # #nn31.benchmark_nlo(Q0=np.sqrt(high2), Q2grid=[low2])
     nn40 = BenchmarkNNPDF40()
-    # nn40.benchmark_nlo(Q2grid=[100])
-    nn40.benchmark_nnlo(Q2grid=[100])
+    # nn40.benchmark_nnlo(Q2grid=[100])
+    nn40.benchmark_nnlo(Q0=np.sqrt(high2), Q2grid=[low2])

benchmarks/asv.conf.json

@@ -0,0 +1,91 @@
+{
+  // The version of the config file format.  Do not change, unless
+  // you know what you are doing.
+  "version": 1,
+
+  // The name of the project being benchmarked
+  "project": "eko",
+
+  // The project's homepage
+  "project_url": "https://n3pdf.github.io/eko/",
+
+  // The URL or local path of the source code repository for the
+  // project being benchmarked
+  "repo": "..",
+
+  // List of branches to benchmark. If not provided, defaults to "master"
+  // (for git) or "tip" (for mercurial).
+  "branches": ["HEAD"],
+
+  // The DVCS being used.  If not set, it will be automatically
+  // determined from "repo" by looking at the protocol in the URL
+  // (if remote), or by looking for special directories, such as
+  // ".git" (if local).
+  "dvcs": "git",
+
+  // The tool to use to create environments.  May be "conda",
+  // "virtualenv" or other value depending on the plugins in use.
+  // If missing or the empty string, the tool will be automatically
+  // determined by looking for tools on the PATH environment
+  // variable.
+  "environment_type": "virtualenv",
+
+  // the base URL to show a commit for the project.
+  "show_commit_url": "https://github.com/N3PDF/eko/commit/",
+
+  // The Pythons you'd like to test against.  If not provided, defaults
+  // to the current version of Python used to run `asv`.
+  "pythons": ["3.8", "3.9", "3.10"],
+
+  // The matrix of dependencies to test.  Each key is the name of a
+  // package (in PyPI) and the values are version numbers.  An empty
+  // list indicates to just test against the default (latest)
+  //version.
+  "matrix": {
+    "poetry": []
+  },
+
+  // The directory (relative to the current directory) that benchmarks are
+  // stored in.  If not provided, defaults to "benchmarks"
+  "benchmark_dir": "performance",
+
+  // The directory (relative to the current directory) to cache the Python
+  // environments in.  If not provided, defaults to "env"
+  // "env_dir": "env",
+
+  // The directory (relative to the current directory) that raw benchmark
+  // results are stored in.  If not provided, defaults to "results".
+  "results_dir": "results",
+
+  // The directory (relative to the current directory) that the html tree
+  // should be written to.  If not provided, defaults to "html".
+  "html_dir": "html",
+
+  // The number of characters to retain in the commit hashes.
+  // "hash_length": 8,
+
+  // `asv` will cache wheels of the recent builds in each
+  // environment, making them faster to install next time.  This is
+  // number of builds to keep, per environment.
+  "build_cache_size": 8,
+
+  // The commits after which the regression search in `asv publish`
+  // should start looking for regressions. Dictionary whose keys are
+  // regexps matching to benchmark names, and values corresponding to
+  // the commit (exclusive) after which to start looking for
+  // regressions.  The default is to start from the first commit
+  // with results. If the commit is `null`, regression detection is
+  // skipped for the matching benchmark.
+  //
+  // "regressions_first_commits": {
+  //    "some_benchmark": "352cdf",  // Consider regressions only after this commit
+  //    "another_benchmark": null,   // Skip regression detection altogether
+  // }
+
+
+  "install_command": ["in-dir={env_dir} python -mpip install {wheel_file}[mark]"],
+
+  "build_command": [
+    "PIP_NO_BUILD_ISOLATION=false python -mpip wheel --no-deps -w {build_cache_dir} {build_dir}"
+  ]
+}

benchmarks/eko/benchmark_ad.py

@@ -4,7 +4,7 @@
 import pytest
 
 import eko.anomalous_dimensions.as2 as ad_as2
-import eko.anomalous_dimensions.harmonics as h
+import eko.harmonics as h
 from eko.constants import CA, CF, TR
 
 
@@ -15,7 +15,7 @@ def benchmark_melling_g3_pegasus():
 
 
 def check_melling_g3_pegasus(N):
-    S1 = h.harmonic_S1(N)
+    S1 = h.S1(N)
     N1 = N + 1.0
     N2 = N + 2.0
     N3 = N + 3.0
@@ -28,7 +28,7 @@ def check_melling_g3_pegasus(N):
     S14 = S13 + 1.0 / N4
     S15 = S14 + 1.0 / N5
     S16 = S15 + 1.0 / N6
-    zeta2 = h.zeta2
+    zeta2 = h.constants.zeta2
 
     SPMOM = (
         1.0000 * (zeta2 - S1 / N) / N
@@ -39,7 +39,7 @@ def check_melling_g3_pegasus(N):
         - 0.3174 * (zeta2 - S15 / N5) / N5
         + 0.0699 * (zeta2 - S16 / N6) / N6
     )
-    np.testing.assert_allclose(h.mellin_g3(N), SPMOM)
+    np.testing.assert_allclose(h.g_functions.mellin_g3(N, S1), SPMOM)
 
 
 @pytest.mark.isolated
@@ -52,12 +52,12 @@ def benchmark_gamma_ns_1_pegasus():
 
 def check_gamma_1_pegasus(N, NF):
     # Test against pegasus implementation
-    ZETA2 = h.zeta2
-    ZETA3 = h.zeta3
+    ZETA2 = h.constants.zeta2
+    ZETA3 = h.constants.zeta3
 
     # N = np.random.rand(1) + np.random.rand(1) * 1j
-    S1 = h.harmonic_S1(N)
-    S2 = h.harmonic_S2(N)
+    S1 = h.S1(N)
+    S2 = h.S2(N)
 
     N1 = N + 1.0
     N2 = N + 2.0
@@ -87,16 +87,18 @@ def check_gamma_1_pegasus(N, NF):
     )
     SLC = -5.0 / 8.0 * ZETA3
     SLV = (
-        -ZETA2 / 2.0 * (h.cern_polygamma(N1 / 2, 0) - h.cern_polygamma(N / 2, 0))
+        -ZETA2
+        / 2.0
+        * (h.polygamma.cern_polygamma(N1 / 2, 0) - h.polygamma.cern_polygamma(N / 2, 0))
         + S1 / NS
         + SPMOM
     )
     SSCHLM = SLC - SLV
-    SSTR2M = ZETA2 - h.cern_polygamma(N1 / 2, 1)
-    SSTR3M = 0.5 * h.cern_polygamma(N1 / 2, 2) + ZETA3
+    SSTR2M = ZETA2 - h.polygamma.cern_polygamma(N1 / 2, 1)
+    SSTR3M = 0.5 * h.polygamma.cern_polygamma(N1 / 2, 2) + ZETA3
     SSCHLP = SLC + SLV
-    SSTR2P = ZETA2 - h.cern_polygamma(N2 / 2, 1)
-    SSTR3P = 0.5 * h.cern_polygamma(N2 / 2, 2) + ZETA3
+    SSTR2P = ZETA2 - h.polygamma.cern_polygamma(N2 / 2, 1)
+    SSTR3P = 0.5 * h.polygamma.cern_polygamma(N2 / 2, 2) + ZETA3
 
     PNMA = (
         16.0 * S1 * (2.0 * N + 1.0) / (NS * N1S)
@@ -136,8 +138,9 @@ def check_gamma_1_pegasus(N, NF):
     P1NSP = CF * ((CF - CA / 2.0) * PNPA + CA * PNSB + TR * NF * PNSC)
     P1NSM = CF * ((CF - CA / 2.0) * PNMA + CA * PNSB + TR * NF * PNSC)
 
-    np.testing.assert_allclose(ad_as2.gamma_nsp(N, NF), -P1NSP)
-    np.testing.assert_allclose(ad_as2.gamma_nsm(N, NF), -P1NSM)
+    sx = h.sx(N, 2)
+    np.testing.assert_allclose(ad_as2.gamma_nsp(N, NF, sx), -P1NSP)
+    np.testing.assert_allclose(ad_as2.gamma_nsm(N, NF, sx), -P1NSM)
 
     NS = N * N
     NT = NS * N
@@ -254,7 +257,7 @@ def check_gamma_1_pegasus(N, NF):
     P1Sgq = (CF * CF * PGQA + CF * CA * PGQB + TR * NF * CF * PGQC) * 4.0
     P1Sgg = (CA * CA * PGGA + TR * NF * (CA * PGGB + CF * PGGC)) * 4.0
 
-    gS1 = ad_as2.gamma_singlet(N, NF)
+    gS1 = ad_as2.gamma_singlet(N, NF, sx)
     np.testing.assert_allclose(gS1[0, 0], -P1Sqq)
     np.testing.assert_allclose(gS1[0, 1], -P1Sqg)
     np.testing.assert_allclose(gS1[1, 0], -P1Sgq)

benchmarks/eko/benchmark_evol_to_unity.py

@@ -49,6 +49,7 @@ class BenchmarkBackwardForward:
         "ev_op_max_order": 1,
         "ev_op_iterations": 1,
         "backward_inversion": "exact",
+        "n_integration_cores": 1,
     }
     g = OperatorGrid.from_dict(
         theory_card,
@@ -99,12 +100,20 @@ def test_operator_grid(
             np.zeros(dim),
             atol=7e-4,
         )
-        # gg
-        np.testing.assert_allclose(
+        # gg, check last two rows separately
+        gg_id = (
             o_back.op_members[(21, 21)].value @ o.op_members[(21, 21)].value
-            + o_back.op_members[(21, 100)].value @ o.op_members[(100, 21)].value,
-            np.eye(dim[0]),
-            atol=9e-2,
+            + o_back.op_members[(21, 100)].value @ o.op_members[(100, 21)].value
+        )
+        np.testing.assert_allclose(
+            gg_id[:-2],
+            np.eye(dim[0])[:-2],
+            atol=3e-2,
+        )
+        np.testing.assert_allclose(
+            gg_id[-2:],
+            np.eye(dim[0])[-2:],
+            atol=11.2e-2,
         )
         # gq
         np.testing.assert_allclose(
@@ -126,14 +135,14 @@ def test_operator_grid(
     #     ome.compute( q2, L)
     #     ome_back.compute(q2, L)
 
-    #     dim = ome.ome_members[(100, 100)].value.shape
+    #     dim = ome.op_members[(100, 100)].value.shape
     #     ome_tensor = np.zeros((3,3,dim[0],dim[0]))
     #     ome_tensor_back = ome_tensor
     #     idx_dict = dict(zip(["g", "q", "H"],[0,1,2]))
     #     for p1, j in idx_dict.items():
     #         for p2, k in idx_dict.items():
-    #             ome_tensor[j,k] = ome.ome_members[f"S_{p1}{p2}"].value
-    #             ome_tensor_back[j,k] = ome_back.ome_members[f"S_{p1}{p2}"].value
+    #             ome_tensor[j,k] = ome.op_members[f"S_{p1}{p2}"].value
+    #             ome_tensor_back[j,k] = ome_back.op_members[f"S_{p1}{p2}"].value
 
     #     ome_product = np.einsum("abjk,bckl -> acjl", ome_tensor_back, ome_tensor)
     #     for j, line in enumerate(ome_product):

benchmarks/lha_paper_bench.py

@@ -222,9 +222,9 @@ def benchmark_sv(self, pto):
 
     # Benchmark to LHA
     # obj = BenchmarkVFNS()
-    # obj = BenchmarkFFNS()
+    obj = BenchmarkFFNS()
 
-    # obj.benchmark_plain(2)
+    # obj.benchmark_plain(0)
     # obj.benchmark_sv(2)
 
     # # VFNS benchmarks with LHA settings

benchmarks/performance/harmonics.py

@@ -0,0 +1,29 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+
+from eko import anomalous_dimensions as ad
+from eko.mellin import Path
+
+NF = 5
+
+
+class TimeSuite:
+    def setup(self):
+        ts = np.linspace(0.5, 1.0 - 1e-5, 100)
+        logx = 0.8
+        axis_offset = True
+        self.ns = []
+        for t in ts:
+            self.ns.append(Path(t, logx, axis_offset).n)
+
+    def time_as1_sing(self):
+        for n in self.ns:
+            ad.gamma_singlet(0, n, NF)
+
+    def time_as2_sing(self):
+        for n in self.ns:
+            ad.gamma_singlet(1, n, NF)
+
+    def time_as3_sing(self):
+        for n in self.ns:
+            ad.gamma_singlet(2, n, NF)

doc/source/code/IO.rst

@@ -12,29 +12,6 @@ The input is split into two runcard
   and environment
 - a operator runcard that defines the numerical setup and the requested operators
 
-.. list-table:: theory input runcard
-  :header-rows: 1
-
-  * - Name
-    - Type
-    - default
-    - description
-  * - ``PTO``
-    - :py:obj:`int`
-    - [required]
-    - order of perturbation theory: ``0`` = LO, ...
-  * - ``alphas``
-    - :py:obj:`float`
-    - [required]
-    - reference value of the strong coupling :math:`\alpha_s(\mu_0^2)` (Note that we have to use
-      :math:`\alpha_s(\mu_0^2)` here, instead of :math:`a_s(\mu_0^2)` for legacy reasons)
-
-.. include:: IO-tabs/ThresholdConfig.rst
-
-.. include:: IO-tabs/InterpolatorDispatcher.rst
-
-- :py:obj:`list(float)` ``Q2grid`` the target grid in :math:`Q^2`; **Required**
-
 
 Output
 ------