Merge pull request #312 from NNPDF/fix_matching_inversion_2

Fix `nf` in inverse matching

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]
+        )

src/eko/evolution_operator/operator_matrix_element.py

@@ -185,8 +185,12 @@ class OperatorMatrixElement(Operator):
         configuration
     managers : dict
         managers
-    segment: Segment
-        path segment
+    nf: int
+        number of active flavor below threshold
+    q2: float
+        squared matching scale
+    is_backward: bool
+        True for backward matching
     L: float
         :math:`\ln(\mu_F^2 / m_h^2)`
     is_msbar: bool

src/eko/runner/parts.py

@@ -152,9 +152,8 @@ def match(eko: EKO, recipe: Matching) -> Operator:
     op.compute()
 
     binfo = blowup_info(eko)
-    nf_match = op.nf - 1 if recipe.inverse else op.nf
     res, err = matching_condition.MatchingCondition.split_ad_to_evol_map(
-        op.op_members, nf_match, recipe.scale, **binfo
+        op.op_members, op.nf, recipe.scale, **binfo
     ).to_flavor_basis_tensor(qed=binfo["qed"])
 
     return Operator(res, err)