From 054e4b2a9c9d9b9ceeaee8aea42963aa5b5c6242 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 31 Jul 2023 12:47:10 +0200 Subject: [PATCH 01/13] Drop max_num_flavs --- benchmarks/eko/benchmark_evol_to_unity.py | 1 - benchmarks/eko/benchmark_strong_coupling.py | 63 ++++++++++--------- benchmarks/ekobox/benchmark_evol_pdf.py | 1 - doc/source/overview/tutorials/alpha_s.ipynb | 7 +-- doc/source/overview/tutorials/pdf.ipynb | 1 - extras/n3lo_bench/splitting_function_utils.py | 1 - src/eko/io/runcards.py | 1 - src/eko/quantities/couplings.py | 1 - src/ekobox/cards.py | 1 - tests/eko/kernels/test_kernels_QEDns.py | 1 - tests/eko/quantities/test_couplings.py | 2 +- tests/eko/test_couplings.py | 5 -- 12 files changed, 37 insertions(+), 48 deletions(-) diff --git a/benchmarks/eko/benchmark_evol_to_unity.py b/benchmarks/eko/benchmark_evol_to_unity.py index 41bc4b4d0..893de5de3 100644 --- a/benchmarks/eko/benchmark_evol_to_unity.py +++ b/benchmarks/eko/benchmark_evol_to_unity.py @@ -20,7 +20,6 @@ def update_cards(theory: TheoryCard, operator: OperatorCard): alphas=0.35, alphaem=0.007496, scale=float(np.sqrt(2)), - max_num_flavs=6, num_flavs_ref=None, ) theory.heavy.num_flavs_init = 4 diff --git a/benchmarks/eko/benchmark_strong_coupling.py b/benchmarks/eko/benchmark_strong_coupling.py index d13a264bd..80f65cf72 100644 --- a/benchmarks/eko/benchmark_strong_coupling.py +++ b/benchmarks/eko/benchmark_strong_coupling.py @@ -8,7 +8,7 @@ from eko.couplings import Couplings from eko.io.runcards import TheoryCard from eko.quantities.couplings import CouplingEvolutionMethod, CouplingsInfo -from eko.quantities.heavy_quarks import QuarkMassScheme +from eko.quantities.heavy_quarks import FlavorsNumber, QuarkMassScheme # try to load LHAPDF - if not available, we'll use the cached values try: @@ -35,16 +35,12 @@ use_PEGASUS = False -def ref_couplings( - ref_values, - ref_scale: float, -) -> CouplingsInfo: +def ref_couplings(ref_values, ref_scale: float, ref_nf: FlavorsNumber) -> CouplingsInfo: return CouplingsInfo( alphas=ref_values[0], alphaem=ref_values[1], scale=ref_scale, - max_num_flavs=6, - num_flavs_ref=None, + num_flavs_ref=ref_nf, ) @@ -54,11 +50,11 @@ def test_a_s(self): """Tests the value of alpha_s (for now only at LO) for a given set of parameters """ - # TODO @JCM: we need a source for this! + # source: JCM known_val = 0.0091807954 ref_mu2 = 90 ask_q2 = 125 - ref = ref_couplings([0.1181, 0.007496], np.sqrt(ref_mu2)) + ref = ref_couplings([0.1181, 0.007496], np.sqrt(ref_mu2), 5) as_FFNS_LO = Couplings( couplings=ref, order=(1, 0), @@ -80,7 +76,7 @@ def benchmark_LHA_paper(self): # LO - FFNS # note that the LO-FFNS value reported in :cite:`Giele:2002hx` # was corrected in :cite:`Dittmar:2005ed` - coupling_ref = ref_couplings([0.35, 0.007496], np.sqrt(2)) + coupling_ref = ref_couplings([0.35, 0.007496], np.sqrt(2), 4) as_FFNS_LO = Couplings( couplings=coupling_ref, order=(1, 0), @@ -141,7 +137,7 @@ def benchmark_APFEL_ffns(self): threshold_holder = matchings.Atlas.ffns(nf, 0.0) for order in [1, 2, 3]: as_FFNS = Couplings( - couplings=ref_couplings(coupling_ref, scale_ref), + couplings=ref_couplings(coupling_ref, scale_ref, nf), order=(order, 0), method=CouplingEvolutionMethod.EXPANDED, masses=threshold_holder.walls[1:-1], @@ -212,8 +208,7 @@ def benchmark_pegasus_ffns(self): } # collect my values threshold_holder = matchings.Atlas.ffns(nf, 0.0) - couplings = ref_couplings(coupling_ref, scale_ref) - couplings.max_num_flavs = 4 + couplings = ref_couplings(coupling_ref, scale_ref, nf) for order in [1, 2, 3, 4]: as_FFNS = Couplings( couplings=couplings, @@ -259,6 +254,7 @@ def benchmark_APFEL_vfns(self): Q2s = [1, 2**2, 3**2, 90**2, 100**2] coupling_ref = np.array([0.118, 0.007496]) scale_ref = 91.0 + nf_ref = 5 threshold_list = np.power([2, 4, 175], 2) apfel_vals_dict = { 1: np.array( @@ -292,7 +288,7 @@ def benchmark_APFEL_vfns(self): # collect my values for order in [1, 2, 3]: as_VFNS = Couplings( - couplings=ref_couplings(coupling_ref, scale_ref), + couplings=ref_couplings(coupling_ref, scale_ref, nf_ref), order=(order, 0), method=CouplingEvolutionMethod.EXPANDED, masses=threshold_list, @@ -339,8 +335,9 @@ def benchmark_APFEL_vfns_fact_to_ren(self): ] coupling_ref = np.array([0.118, 0.007496]) scale_ref = 91.0 + nf_ref = 5 fact_to_ren_lin_list = [0.567, 2.34] - threshold_list = np.power([2, 2 * 4, 2 * 92], 2) + masses_list = np.power([2, 2 * 4, 2 * 92], 2) apfel_vals_dict_list = [ { 1: np.array( @@ -437,10 +434,10 @@ def benchmark_APFEL_vfns_fact_to_ren(self): # collect my values for order in [1, 2, 3]: as_VFNS = Couplings( - couplings=ref_couplings(coupling_ref, scale_ref), + couplings=ref_couplings(coupling_ref, scale_ref, nf_ref), order=(order, 0), method=CouplingEvolutionMethod.EXACT, - masses=threshold_list.tolist(), + masses=masses_list.tolist(), hqm_scheme=QuarkMassScheme.POLE, thresholds_ratios=np.array([1.0, 1.0, 1.0]) / fact_to_ren_lin**2, ) @@ -457,7 +454,7 @@ def benchmark_APFEL_vfns_fact_to_ren(self): apfel.SetAlphaEvolution("exact") apfel.SetAlphaQCDRef(coupling_ref[0], scale_ref) apfel.SetVFNS() - apfel.SetPoleMasses(*np.sqrt(threshold_list)) + apfel.SetPoleMasses(*np.sqrt(masses_list)) apfel.SetRenFacRatio(1.0 / fact_to_ren_lin) apfel.InitializeAPFEL() # collect a_s @@ -468,12 +465,18 @@ def benchmark_APFEL_vfns_fact_to_ren(self): ) 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=2.5e-5) + np.testing.assert_allclose( + apfel_vals, + np.array(my_vals), + rtol=2.5e-5, + err_msg=f"{order=},{fact_to_ren_lin=}", + ) def benchmark_APFEL_vfns_threshold(self): Q2s = np.power([30, 96, 150], 2) coupling_ref = np.array([0.118, 0.007496]) scale_ref = 91.0 + nf_ref = 4 threshold_list = np.power([30, 95, 240], 2) thresholds_ratios = [2.34**2, 1.0**2, 0.5**2] apfel_vals_dict = { @@ -487,7 +490,7 @@ def benchmark_APFEL_vfns_threshold(self): # collect my values for order in [2, 3]: as_VFNS = Couplings( - couplings=ref_couplings(coupling_ref, scale_ref), + couplings=ref_couplings(coupling_ref, scale_ref, nf_ref), order=(order, 0), method=CouplingEvolutionMethod.EXPANDED, masses=threshold_list.tolist(), @@ -496,7 +499,6 @@ def benchmark_APFEL_vfns_threshold(self): ) my_vals = [] for Q2 in Q2s: - print(Q2) my_vals.append(as_VFNS.a(Q2)[0]) # get APFEL numbers - if available else use cache apfel_vals = apfel_vals_dict[order] @@ -516,7 +518,6 @@ def benchmark_APFEL_vfns_threshold(self): 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)) @@ -525,6 +526,7 @@ def benchmark_APFEL_vfns_msbar(self): Q2s = np.power([3, 96, 150], 2) coupling_ref = np.array([0.118, 0.007496]) scale_ref = 91.0 + nf_ref = 5 Q2m = np.power([2.0, 4.0, 175], 2) m2s = np.power((1.4, 4.0, 175), 2) apfel_vals_dict = { @@ -538,7 +540,7 @@ def benchmark_APFEL_vfns_msbar(self): # collect my values for order in [2, 3]: as_VFNS = Couplings( - couplings=ref_couplings(coupling_ref, scale_ref), + couplings=ref_couplings(coupling_ref, scale_ref, nf_ref), order=(order, 0), method=CouplingEvolutionMethod.EXPANDED, masses=m2s.tolist(), @@ -586,7 +588,7 @@ def benchmark_lhapdf_ffns_lo(self): # collect my values threshold_holder = matchings.Atlas.ffns(nf, 0.0) as_FFNS_LO = Couplings( - couplings=ref_couplings(coupling_ref, scale_ref), + couplings=ref_couplings(coupling_ref, scale_ref, nf), order=(1, 0), method=CouplingEvolutionMethod.EXACT, masses=threshold_holder.walls[1:-1], @@ -629,8 +631,9 @@ def benchmark_apfel_exact(self): Q2s = [1e1, 1e2, 1e3, 1e4] coupling_ref = np.array([0.118, 0.007496]) scale_ref = 90 + nf = 3 # collect my values - threshold_holder = matchings.Atlas.ffns(3, 0.0) + threshold_holder = matchings.Atlas.ffns(nf, 0.0) # LHAPDF cache apfel_vals_dict = { 1: np.array( @@ -660,7 +663,7 @@ def benchmark_apfel_exact(self): } for order in range(1, 3 + 1): sc = Couplings( - couplings=ref_couplings(coupling_ref, scale_ref), + couplings=ref_couplings(coupling_ref, scale_ref, nf), order=(order, 0), method=CouplingEvolutionMethod.EXACT, masses=threshold_holder.walls[1:-1], @@ -696,8 +699,9 @@ def benchmark_lhapdf_exact(self): Q2s = [1e1, 1e2, 1e3, 1e4] coupling_ref = np.array([0.118, 0.007496]) scale_ref = 90 + nf = 3 # collect my values - threshold_holder = matchings.Atlas.ffns(3, 0.0) + threshold_holder = matchings.Atlas.ffns(nf, 0.0) # LHAPDF cache lhapdf_vals_dict = { 1: np.array( @@ -735,7 +739,7 @@ def benchmark_lhapdf_exact(self): } for order in range(1, 4 + 1): sc = Couplings( - couplings=ref_couplings(coupling_ref, scale_ref), + couplings=ref_couplings(coupling_ref, scale_ref, nf), order=(order, 0), method=CouplingEvolutionMethod.EXACT, masses=threshold_holder.walls[1:-1], @@ -766,6 +770,7 @@ def benchmark_lhapdf_zmvfns_lo(self): Q2s = [1, 1e1, 1e2, 1e3, 1e4] coupling_ref = np.array([0.118, 0.007496]) scale_ref = 900 + nf_ref = 5 m2c = 2.0 m2b = 25.0 m2t = 1500.0 @@ -785,7 +790,7 @@ def benchmark_lhapdf_zmvfns_lo(self): # collect my values as_VFNS_LO = Couplings( - couplings=ref_couplings(coupling_ref, np.sqrt(scale_ref)), + couplings=ref_couplings(coupling_ref, np.sqrt(scale_ref), nf_ref), order=(1, 0), method=CouplingEvolutionMethod.EXACT, masses=threshold_list, diff --git a/benchmarks/ekobox/benchmark_evol_pdf.py b/benchmarks/ekobox/benchmark_evol_pdf.py index 3a9dc0f4a..2a9aa8ea9 100644 --- a/benchmarks/ekobox/benchmark_evol_pdf.py +++ b/benchmarks/ekobox/benchmark_evol_pdf.py @@ -23,7 +23,6 @@ def benchmark_evolve_single_member( theory.couplings.alphas = 0.118000 theory.couplings.scale = 91.1876 theory.couplings.alphaem = 0.007496 - theory.couplings.max_num_flavs = 3 theory.heavy.num_flavs_max_pdf = 3 theory.heavy.masses.c.value = 1.3 theory.heavy.masses.b.value = 4.75 diff --git a/doc/source/overview/tutorials/alpha_s.ipynb b/doc/source/overview/tutorials/alpha_s.ipynb index a9d75c0ca..d8954900d 100644 --- a/doc/source/overview/tutorials/alpha_s.ipynb +++ b/doc/source/overview/tutorials/alpha_s.ipynb @@ -40,7 +40,7 @@ "\n", "# set the (alpha_s, alpha_em) reference values\n", "couplings_ref = CouplingsInfo(\n", - " alphas=0.118, alphaem=0.007496252, scale=91.0, num_flavs_ref=None, max_num_flavs=5\n", + " alphas=0.118, alphaem=0.007496252, scale=91.0, num_flavs_ref=5\n", ")\n", "\n", "# set heavy quark masses and their threshold ratios\n", @@ -105,13 +105,10 @@ } ], "metadata": { - "interpreter": { - "hash": "0a84ba3ac8703c04e87bc503a7d00188dfd591ad56130da93c406115a1e4a408" - }, "kernelspec": { "display_name": "eko-KkPVjVhh-py3.10", "language": "python", - "name": "eko-kkpvjvhh-py3.10" + "name": "python3" }, "language_info": { "codemirror_mode": { diff --git a/doc/source/overview/tutorials/pdf.ipynb b/doc/source/overview/tutorials/pdf.ipynb index 4607ef891..9da9ead26 100644 --- a/doc/source/overview/tutorials/pdf.ipynb +++ b/doc/source/overview/tutorials/pdf.ipynb @@ -371,7 +371,6 @@ "th_card.couplings.alphas = 0.130000 # reference value of alpha_s\n", "th_card.couplings.scale = 91.1876 # the reference scale at which alpha_s is provided\n", "th_card.couplings.num_flavs_ref = 5 # the number of flavors active at the alpha_s reference scale\n", - "th_card.couplings.max_num_flavs = 5 # the maximum number of flavors active in the alpha_s evolution\n", "th_card.couplings.num_flavs_init = 3 # the number of flavors active at the reference scale\n", "th_card.num_flavs_max_pdf = 5 # the maximum number of flavors active in the pdf evolution." ] diff --git a/extras/n3lo_bench/splitting_function_utils.py b/extras/n3lo_bench/splitting_function_utils.py index 6917e6bd8..db06f5c72 100644 --- a/extras/n3lo_bench/splitting_function_utils.py +++ b/extras/n3lo_bench/splitting_function_utils.py @@ -104,7 +104,6 @@ def compute_a_s(q2=None, xif2=1.0, nf=None, order=(4, 0)): alphas=0.1181, alphaem=0.007496, scale=91.00, - max_num_flavs=6, num_flavs_ref=5, ) sc = Couplings( diff --git a/src/eko/io/runcards.py b/src/eko/io/runcards.py index 41aef9445..6f1b585b0 100644 --- a/src/eko/io/runcards.py +++ b/src/eko/io/runcards.py @@ -185,7 +185,6 @@ def new_theory(self): em_running=em_running, scale=old["Qref"], num_flavs_ref=old["nfref"], - max_num_flavs=old["MaxNfAs"], ) new["heavy"] = { "num_flavs_init": nf_default(old["Q0"] ** 2.0, default_atlas(ms, ks)) diff --git a/src/eko/quantities/couplings.py b/src/eko/quantities/couplings.py index 01b11d138..f80420a9c 100644 --- a/src/eko/quantities/couplings.py +++ b/src/eko/quantities/couplings.py @@ -20,7 +20,6 @@ class CouplingsInfo(DictLike): alphas: Coupling alphaem: Coupling scale: LinearScale - max_num_flavs: FlavorsNumber num_flavs_ref: FlavorsNumber r"""Number of active flavors at strong coupling reference scale. diff --git a/src/ekobox/cards.py b/src/ekobox/cards.py index 474e7d931..e68a8ce63 100644 --- a/src/ekobox/cards.py +++ b/src/ekobox/cards.py @@ -15,7 +15,6 @@ alphaem=0.007496252, scale=91.2, num_flavs_ref=5, - max_num_flavs=6, ), heavy=dict( num_flavs_init=4, diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py index e84956084..b85a3f013 100644 --- a/tests/eko/kernels/test_kernels_QEDns.py +++ b/tests/eko/kernels/test_kernels_QEDns.py @@ -121,7 +121,6 @@ def test_zero_true_gamma(): alphaem=alpharef[1], scale=muref, num_flavs_ref=5, - max_num_flavs=6, ) ) evmod = CouplingEvolutionMethod.EXACT diff --git a/tests/eko/quantities/test_couplings.py b/tests/eko/quantities/test_couplings.py index 00ccae9d1..83bef709c 100644 --- a/tests/eko/quantities/test_couplings.py +++ b/tests/eko/quantities/test_couplings.py @@ -3,7 +3,7 @@ def test_couplings_ref(): scale = 90.0 - d = dict(alphas=0.1, alphaem=0.01, scale=scale, max_num_flavs=6, num_flavs_ref=5) + d = dict(alphas=0.1, alphaem=0.01, scale=scale, num_flavs_ref=5) couplings = CouplingsInfo.from_dict(d) assert couplings.scale == scale assert not couplings.em_running diff --git a/tests/eko/test_couplings.py b/tests/eko/test_couplings.py index 646a56448..10535cbcb 100644 --- a/tests/eko/test_couplings.py +++ b/tests/eko/test_couplings.py @@ -55,7 +55,6 @@ def test_init(self): alphaem=alpharef[1], scale=muref, num_flavs_ref=5, - max_num_flavs=6, ) ) order = (1, 0) @@ -162,7 +161,6 @@ def test_ref(self): alphaem=alpharef[1], scale=muref, num_flavs_ref=nfref, - max_num_flavs=6, ) ) for order_s in [1, 2, 3, 4]: @@ -197,7 +195,6 @@ def test_ref_copy_e841b0dfdee2f31d9ccc1ecee4d9d1a6f6624313(self): alphaem=alpharef[1], scale=muref, num_flavs_ref=3, # reference nf is needed to force the matching - max_num_flavs=6, ) ) sc = Couplings( @@ -238,7 +235,6 @@ def test_exact(self): alphaem=alpharef[1], scale=muref, num_flavs_ref=nfref, - max_num_flavs=6, em_running=em_running, ) ) @@ -301,7 +297,6 @@ def benchmark_expanded_n3lo(self): alphaem=alpharef[1], scale=muref, num_flavs_ref=5, - max_num_flavs=6, ) m2c = 2 m2b = 25 From cfc7ee53bd807c5377263fd1024e313d4d5f8315 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 31 Jul 2023 12:58:29 +0200 Subject: [PATCH 02/13] Drop num_flavs_max_pdf --- benchmarks/ekobox/benchmark_evol_pdf.py | 1 - doc/source/overview/tutorials/pdf.ipynb | 3 +-- src/eko/io/legacy.py | 1 - src/eko/io/runcards.py | 1 - src/eko/quantities/heavy_quarks.py | 2 -- src/ekobox/cards.py | 1 - tests/eko/test_quantities.py | 1 - 7 files changed, 1 insertion(+), 9 deletions(-) diff --git a/benchmarks/ekobox/benchmark_evol_pdf.py b/benchmarks/ekobox/benchmark_evol_pdf.py index 2a9aa8ea9..c85083604 100644 --- a/benchmarks/ekobox/benchmark_evol_pdf.py +++ b/benchmarks/ekobox/benchmark_evol_pdf.py @@ -23,7 +23,6 @@ def benchmark_evolve_single_member( theory.couplings.alphas = 0.118000 theory.couplings.scale = 91.1876 theory.couplings.alphaem = 0.007496 - theory.heavy.num_flavs_max_pdf = 3 theory.heavy.masses.c.value = 1.3 theory.heavy.masses.b.value = 4.75 theory.heavy.masses.t.value = 172 diff --git a/doc/source/overview/tutorials/pdf.ipynb b/doc/source/overview/tutorials/pdf.ipynb index 9da9ead26..0719d9d63 100644 --- a/doc/source/overview/tutorials/pdf.ipynb +++ b/doc/source/overview/tutorials/pdf.ipynb @@ -371,8 +371,7 @@ "th_card.couplings.alphas = 0.130000 # reference value of alpha_s\n", "th_card.couplings.scale = 91.1876 # the reference scale at which alpha_s is provided\n", "th_card.couplings.num_flavs_ref = 5 # the number of flavors active at the alpha_s reference scale\n", - "th_card.couplings.num_flavs_init = 3 # the number of flavors active at the reference scale\n", - "th_card.num_flavs_max_pdf = 5 # the maximum number of flavors active in the pdf evolution." + "th_card.heavy.num_flavs_init = 3 # the number of flavors active at the reference scale" ] }, { diff --git a/src/eko/io/legacy.py b/src/eko/io/legacy.py index d580f7373..cb7b24644 100644 --- a/src/eko/io/legacy.py +++ b/src/eko/io/legacy.py @@ -100,7 +100,6 @@ def from_old(cls, old: RawCard): """Load from old metadata.""" heavy = HeavyInfo( num_flavs_init=4, - num_flavs_max_pdf=5, intrinsic_flavors=[], masses=HeavyQuarkMasses( [ diff --git a/src/eko/io/runcards.py b/src/eko/io/runcards.py index 6f1b585b0..c66e48dc4 100644 --- a/src/eko/io/runcards.py +++ b/src/eko/io/runcards.py @@ -190,7 +190,6 @@ def new_theory(self): "num_flavs_init": nf_default(old["Q0"] ** 2.0, default_atlas(ms, ks)) if old["nf0"] is None else old["nf0"], - "num_flavs_max_pdf": old["MaxNfPdf"], "matching_ratios": self.heavies("k%sThr", old), "masses_scheme": old["HQ"], } diff --git a/src/eko/quantities/heavy_quarks.py b/src/eko/quantities/heavy_quarks.py index df0285649..ab0ece64a 100644 --- a/src/eko/quantities/heavy_quarks.py +++ b/src/eko/quantities/heavy_quarks.py @@ -96,8 +96,6 @@ class HeavyInfo(DictLike): I.e. :math:`n_{f,\text{ref}}(\mu^2_0)`, formerly called ``nf0``. """ - num_flavs_max_pdf: FlavorsNumber - """Maximum number of quark PDFs.""" intrinsic_flavors: IntrinsicFlavors """List of intrinsic quark PDFs.""" masses: HeavyQuarkMasses diff --git a/src/ekobox/cards.py b/src/ekobox/cards.py index e68a8ce63..a69b32ad6 100644 --- a/src/ekobox/cards.py +++ b/src/ekobox/cards.py @@ -18,7 +18,6 @@ ), heavy=dict( num_flavs_init=4, - num_flavs_max_pdf=6, intrinsic_flavors=[4], masses=[ReferenceRunning([mq, nan]) for mq in (2.0, 4.5, 173.07)], masses_scheme="POLE", diff --git a/tests/eko/test_quantities.py b/tests/eko/test_quantities.py index 0c2628b30..8fd6c298f 100644 --- a/tests/eko/test_quantities.py +++ b/tests/eko/test_quantities.py @@ -34,7 +34,6 @@ def test_HeavyQuarks(): def test_HeavyInfo(): i = hq.HeavyInfo( num_flavs_init=4, - num_flavs_max_pdf=6, intrinsic_flavors=[4, 5], masses=hq.HeavyQuarkMasses( [ From 8b3899d02dd80e6ea4dc31a0c517ac7785bcc3d8 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 31 Jul 2023 19:13:24 +0200 Subject: [PATCH 03/13] Start dropping intrinsic --- benchmarks/eko/benchmark_evol_to_unity.py | 1 - src/eko/evolution_operator/grid.py | 8 +- .../evolution_operator/matching_condition.py | 48 +++++----- .../operator_matrix_element.py | 8 +- src/eko/evolution_operator/physical.py | 23 ++--- src/eko/io/legacy.py | 1 - src/eko/io/runcards.py | 5 - src/eko/io/types.py | 1 - src/eko/quantities/heavy_quarks.py | 10 +- src/eko/runner/legacy.py | 2 - src/eko/runner/managed.py | 2 - src/eko/runner/parts.py | 40 ++------ src/ekobox/cards.py | 1 - .../test_matching_condition.py | 91 +++++-------------- tests/eko/evolution_operator/test_physical.py | 36 +++----- tests/eko/io/test_runcards.py | 1 - tests/eko/test_quantities.py | 1 - 17 files changed, 80 insertions(+), 199 deletions(-) diff --git a/benchmarks/eko/benchmark_evol_to_unity.py b/benchmarks/eko/benchmark_evol_to_unity.py index 893de5de3..d269db202 100644 --- a/benchmarks/eko/benchmark_evol_to_unity.py +++ b/benchmarks/eko/benchmark_evol_to_unity.py @@ -23,7 +23,6 @@ def update_cards(theory: TheoryCard, operator: OperatorCard): num_flavs_ref=None, ) theory.heavy.num_flavs_init = 4 - theory.heavy.intrinsic_flavors = [4, 5] theory.heavy.masses.c.value = 1.0 theory.heavy.masses.b.value = 4.75 theory.heavy.masses.t.value = 173.0 diff --git a/src/eko/evolution_operator/grid.py b/src/eko/evolution_operator/grid.py index 1180697a6..382ece751 100644 --- a/src/eko/evolution_operator/grid.py +++ b/src/eko/evolution_operator/grid.py @@ -60,7 +60,6 @@ def __init__( masses: List[float], mass_scheme, thresholds_ratios: List[float], - intrinsic_flavors: list, xif: float, n3lo_ad_variation: tuple, configs: Configs, @@ -72,7 +71,6 @@ def __init__( # check config: Dict[str, Any] = {} config["order"] = order - config["intrinsic_range"] = intrinsic_flavors config["xif2"] = xif**2 config["HQ"] = mass_scheme config["ModSV"] = configs.scvar_method @@ -171,16 +169,15 @@ def generate(self, q2: EPoint) -> OpDict: operator.compute() is_downward = is_downward_path(path) - intrinsic_range = [4, 5, 6] if is_downward else self.config["intrinsic_range"] qed = self.config["order"][1] > 0 final_op = physical.PhysicalOperator.ad_to_evol_map( - operator.op_members, operator.nf, operator.q2_to, intrinsic_range, qed + operator.op_members, operator.nf, operator.q2_to, qed ) # and multiply the lower ones from the right for op in reversed(list(thr_ops)): phys_op = physical.PhysicalOperator.ad_to_evol_map( - op.op_members, op.nf, op.q2_to, intrinsic_range, qed + op.op_members, op.nf, op.q2_to, qed ) # join with the basis rotation, since matching requires c+ (or likewise) @@ -189,7 +186,6 @@ def generate(self, q2: EPoint) -> OpDict: self._matching_operators[op.q2_to], nf_match, op.q2_to, - intrinsic_range=intrinsic_range, qed=qed, ) if is_downward: diff --git a/src/eko/evolution_operator/matching_condition.py b/src/eko/evolution_operator/matching_condition.py index 059617fbf..e84ef4d2e 100644 --- a/src/eko/evolution_operator/matching_condition.py +++ b/src/eko/evolution_operator/matching_condition.py @@ -20,7 +20,6 @@ def split_ad_to_evol_map( op_members, nf, q2_thr, - intrinsic_range, qed=False, ): """ @@ -35,8 +34,6 @@ def split_ad_to_evol_map( number of active flavors *below* the threshold q2_thr: float threshold value - intrinsic_range : list - list of intrinsic quark pids qed : bool activate qed """ @@ -78,27 +75,26 @@ def split_ad_to_evol_map( ) # intrinsic matching - if len(intrinsic_range) != 0: - op_id = member.OpMember.id_like(op_members[(200, 200)]) - for intr_fl in intrinsic_range: - ihq = br.quark_names[intr_fl - 1] # find name - if intr_fl > nf + 1: - # keep the higher quarks as they are - m[f"{ihq}+.{ihq}+"] = op_id.copy() - m[f"{ihq}-.{ihq}-"] = op_id.copy() - elif intr_fl == nf + 1: - # match the missing contribution from h+ and h- - m.update( - { - f"{ihq}+.{ihq}+": op_members[ - (br.matching_hplus_pid, br.matching_hplus_pid) - ], - f"S.{ihq}+": op_members[(100, br.matching_hplus_pid)], - f"g.{ihq}+": op_members[(21, br.matching_hplus_pid)], - f"{ihq}-.{ihq}-": op_members[ - (br.matching_hminus_pid, br.matching_hminus_pid) - ], - # f"V.{ihq}-": op_members[(200, br.matching_hminus_pid)], - } - ) + op_id = member.OpMember.id_like(op_members[(200, 200)]) + for intr_fl in [4, 5, 6]: + ihq = br.quark_names[intr_fl - 1] # find name + if intr_fl > nf + 1: + # keep the higher quarks as they are + m[f"{ihq}+.{ihq}+"] = op_id.copy() + m[f"{ihq}-.{ihq}-"] = op_id.copy() + elif intr_fl == nf + 1: + # match the missing contribution from h+ and h- + m.update( + { + f"{ihq}+.{ihq}+": op_members[ + (br.matching_hplus_pid, br.matching_hplus_pid) + ], + f"S.{ihq}+": op_members[(100, br.matching_hplus_pid)], + f"g.{ihq}+": op_members[(21, br.matching_hplus_pid)], + f"{ihq}-.{ihq}-": op_members[ + (br.matching_hminus_pid, br.matching_hminus_pid) + ], + # f"V.{ihq}-": op_members[(200, br.matching_hminus_pid)], + } + ) return cls.promote_names(m, q2_thr) diff --git a/src/eko/evolution_operator/operator_matrix_element.py b/src/eko/evolution_operator/operator_matrix_element.py index 8966a62b9..b45e06c81 100644 --- a/src/eko/evolution_operator/operator_matrix_element.py +++ b/src/eko/evolution_operator/operator_matrix_element.py @@ -214,10 +214,6 @@ class OperatorMatrixElement(Operator): def __init__(self, config, managers, nf, q2, is_backward, L, is_msbar): super().__init__(config, managers, Segment(q2, q2, nf)) self.backward_method = config["backward_inversion"] if is_backward else None - if is_backward: - self.is_intrinsic = True - else: - self.is_intrinsic = bool(len(config["intrinsic_range"]) != 0) self.L = L self.is_msbar = is_msbar # Note for the moment only QCD matching is implemented @@ -238,7 +234,7 @@ def labels(self): logger.warning("%s: skipping non-singlet sector", self.log_label) else: labels.append((200, 200)) - if self.is_intrinsic or self.backward_method is not None: + if self.backward_method is not None: # intrinsic labels, which are not zero at NLO labels.append((br.matching_hminus_pid, br.matching_hminus_pid)) # These contributions are always 0 for the moment @@ -254,7 +250,7 @@ def labels(self): (br.matching_hplus_pid, 100), ] ) - if self.is_intrinsic or self.backward_method is not None: + if self.backward_method is not None: labels.extend( [ (21, br.matching_hplus_pid), diff --git a/src/eko/evolution_operator/physical.py b/src/eko/evolution_operator/physical.py index c963e2a26..6a5138d25 100644 --- a/src/eko/evolution_operator/physical.py +++ b/src/eko/evolution_operator/physical.py @@ -22,7 +22,7 @@ class PhysicalOperator(member.OperatorBase): """ @classmethod - def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed=False): + def ad_to_evol_map(cls, op_members, nf, q2_final, qed=False): """ Obtain map between the 3-dimensional anomalous dimension basis and the 4-dimensional evolution basis. @@ -34,8 +34,6 @@ def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed=False): operator members in anomalous dimension basis nf : int number of active light flavors - intrinsic_range : sequence - intrinsic heavy flavors qed : bool activate qed @@ -93,15 +91,14 @@ def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed=False): m["Tu8.Tu8"] = op_members[(br.non_singlet_pids_map["ns+u"], 0)] m["Vu8.Vu8"] = op_members[(br.non_singlet_pids_map["ns-u"], 0)] # deal with intrinsic heavy quark PDFs - if intrinsic_range is not None: - hqfl = "cbt" - op_id = member.OpMember.id_like(op_members[(21, 21)]) - for intr_fl in intrinsic_range: - if intr_fl <= nf: # light quarks are not intrinsic - continue - hq = hqfl[intr_fl - 4] # find name - # intrinsic means no evolution, i.e. they are evolving with the identity - m[f"{hq}+.{hq}+"] = op_id.copy() - m[f"{hq}-.{hq}-"] = op_id.copy() + hqfl = "cbt" + op_id = member.OpMember.id_like(op_members[(21, 21)]) + for intr_fl in [4, 5, 6]: + if intr_fl <= nf: # light quarks are not intrinsic + continue + hq = hqfl[intr_fl - 4] # find name + # intrinsic means no evolution, i.e. they are evolving with the identity + m[f"{hq}+.{hq}+"] = op_id.copy() + m[f"{hq}-.{hq}-"] = op_id.copy() # map key to MemberName return cls.promote_names(m, q2_final) diff --git a/src/eko/io/legacy.py b/src/eko/io/legacy.py index cb7b24644..8491a6acc 100644 --- a/src/eko/io/legacy.py +++ b/src/eko/io/legacy.py @@ -100,7 +100,6 @@ def from_old(cls, old: RawCard): """Load from old metadata.""" heavy = HeavyInfo( num_flavs_init=4, - intrinsic_flavors=[], masses=HeavyQuarkMasses( [ ReferenceRunning([_MC, np.inf]), diff --git a/src/eko/io/runcards.py b/src/eko/io/runcards.py index c66e48dc4..d5c8c36e8 100644 --- a/src/eko/io/runcards.py +++ b/src/eko/io/runcards.py @@ -193,11 +193,6 @@ def new_theory(self): "matching_ratios": self.heavies("k%sThr", old), "masses_scheme": old["HQ"], } - intrinsic = [] - for idx, q in enumerate(hq.FLAVORS): - if old.get(f"i{q}".upper()) == 1: - intrinsic.append(idx + 4) - new["heavy"]["intrinsic_flavors"] = intrinsic if old["HQ"] == "POLE": new["heavy"]["masses"] = [[m, nan] for m in ms] elif old["HQ"] == "MSBAR": diff --git a/src/eko/io/types.py b/src/eko/io/types.py index 5573bafbb..b92628116 100644 --- a/src/eko/io/types.py +++ b/src/eko/io/types.py @@ -22,7 +22,6 @@ Order = Tuple[int, int] FlavorsNumber = int FlavorIndex = int -IntrinsicFlavors = typing.List[FlavorIndex] N3LOAdVariation = typing.Tuple[int, int, int, int] # Evolution coordinates diff --git a/src/eko/quantities/heavy_quarks.py b/src/eko/quantities/heavy_quarks.py index ab0ece64a..3b081a03d 100644 --- a/src/eko/quantities/heavy_quarks.py +++ b/src/eko/quantities/heavy_quarks.py @@ -6,13 +6,7 @@ import numpy as np from ..io.dictlike import DictLike -from ..io.types import ( - FlavorsNumber, - IntrinsicFlavors, - LinearScale, - ReferenceRunning, - SquaredScale, -) +from ..io.types import FlavorsNumber, LinearScale, ReferenceRunning, SquaredScale FLAVORS = "cbt" @@ -96,8 +90,6 @@ class HeavyInfo(DictLike): I.e. :math:`n_{f,\text{ref}}(\mu^2_0)`, formerly called ``nf0``. """ - intrinsic_flavors: IntrinsicFlavors - """List of intrinsic quark PDFs.""" masses: HeavyQuarkMasses """List of heavy quark masses.""" masses_scheme: QuarkMassScheme diff --git a/src/eko/runner/legacy.py b/src/eko/runner/legacy.py index c964c3415..76d6f1a8a 100644 --- a/src/eko/runner/legacy.py +++ b/src/eko/runner/legacy.py @@ -44,7 +44,6 @@ def __init__( """ new_theory, new_operator = runcards.update(theory_card, operators_card) - new_theory.heavy.intrinsic_flavors = [4, 5, 6] # Store inputs self.path = path @@ -69,7 +68,6 @@ def __init__( masses=masses, mass_scheme=new_theory.heavy.masses_scheme.value, thresholds_ratios=new_theory.heavy.squared_ratios, - intrinsic_flavors=new_theory.heavy.intrinsic_flavors, xif=new_theory.xif, configs=new_operator.configs, debug=new_operator.debug, diff --git a/src/eko/runner/managed.py b/src/eko/runner/managed.py index c353decde..dc253e2ea 100644 --- a/src/eko/runner/managed.py +++ b/src/eko/runner/managed.py @@ -20,8 +20,6 @@ def solve(theory: TheoryCard, operator: OperatorCard, path: Path): """Solve DGLAP equations in terms of evolution kernel operators (EKO).""" - theory.heavy.intrinsic_flavors = [4, 5, 6] - with EKO.create(path) as builder: eko = builder.load_cards(theory, operator).build() # pylint: disable=E1101 diff --git a/src/eko/runner/parts.py b/src/eko/runner/parts.py index 9f97a6de5..22a2ff9ca 100644 --- a/src/eko/runner/parts.py +++ b/src/eko/runner/parts.py @@ -40,32 +40,6 @@ def managers(eko: EKO) -> Managers: ) -def blowup_info(eko: EKO) -> dict: - """Prepare common information to blow up to flavor basis. - - Note - ---- - ``intrinsic_range`` is a fully deprecated feature, here and anywhere else, - since a full range is already always used for backward evolution, and it is - not harmful to use it also for forward. - - Indeed, the only feature of non-intrinsic evolution is to absorb a - non-trivial boundary condition when an intrinsic PDF is defined. - But to achieve this, is sufficient to not specify any intrinsic boundary - condition at all, while if something is there, it is intuitive enough that - it will be consistently evolved. - - Moreover, since two different behavior are applied for the forward and - backward evolution, the intrinsic range is a "non-local" function, since it - does not depend only on the evolution segment, but also on the previous - evolution history (to determine if evolution is backward in flavor, - irrespectively of happening for an increasing or decreasing interval in - scale at fixed flavor). - - """ - return dict(intrinsic_range=[4, 5, 6], qed=eko.theory_card.order[1] > 0) - - def evolve_configs(eko: EKO) -> dict: """Create configs for :class:`Operator`. @@ -99,10 +73,10 @@ def evolve(eko: EKO, recipe: Evolution) -> Operator: ) op.compute() - binfo = blowup_info(eko) + qed = eko.theory_card.order[1] > 0 res, err = physical.PhysicalOperator.ad_to_evol_map( - op.op_members, op.nf, op.q2_to, **binfo - ).to_flavor_basis_tensor(qed=binfo["qed"]) + op.op_members, op.nf, op.q2_to, qed + ).to_flavor_basis_tensor(qed) return Operator(res, err) @@ -120,7 +94,6 @@ def matching_configs(eko: EKO) -> dict: return dict( **evolve_configs(eko), backward_inversion=ocard.configs.inversion_method, - intrinsic_range=tcard.heavy.intrinsic_flavors, ) @@ -148,11 +121,10 @@ def match(eko: EKO, recipe: Matching) -> Operator: eko.theory_card.heavy.masses_scheme is QuarkMassScheme.MSBAR, ) op.compute() - - binfo = blowup_info(eko) + qed = eko.theory_card.order[1] > 0 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 - ).to_flavor_basis_tensor(qed=binfo["qed"]) + op.op_members, nf_match, recipe.scale, qed + ).to_flavor_basis_tensor(qed) return Operator(res, err) diff --git a/src/ekobox/cards.py b/src/ekobox/cards.py index a69b32ad6..482e65faf 100644 --- a/src/ekobox/cards.py +++ b/src/ekobox/cards.py @@ -18,7 +18,6 @@ ), heavy=dict( num_flavs_init=4, - intrinsic_flavors=[4], masses=[ReferenceRunning([mq, nan]) for mq in (2.0, 4.5, 173.07)], masses_scheme="POLE", matching_ratios=[1.0, 1.0, 1.0], diff --git a/tests/eko/evolution_operator/test_matching_condition.py b/tests/eko/evolution_operator/test_matching_condition.py index 8cbd6b537..bb6aeebbc 100644 --- a/tests/eko/evolution_operator/test_matching_condition.py +++ b/tests/eko/evolution_operator/test_matching_condition.py @@ -22,12 +22,6 @@ def mkOME(self): (br.matching_hplus_pid, 21), (200, 200), (br.matching_hminus_pid, 200), - ]: - ome.update({key: mkOM(self.shape)}) - return ome - - def update_intrinsic_OME(self, ome): - for key in [ (br.matching_hplus_pid, br.matching_hplus_pid), (br.matching_hminus_pid, br.matching_hminus_pid), (200, br.matching_hminus_pid), @@ -35,10 +29,11 @@ def update_intrinsic_OME(self, ome): (21, br.matching_hplus_pid), ]: ome.update({key: mkOM(self.shape)}) + return ome def test_split_ad_to_evol_map(self): ome = self.mkOME() - a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, []) + a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1) triv_keys = [ "V.V", "T3.T3", @@ -55,6 +50,14 @@ def test_split_ad_to_evol_map(self): "c+.S", "c+.g", # "c-.V", + "S.c+", + "g.c+", + "c+.c+", + "c-.c-", + "b+.b+", + "b-.b-", + "t+.t+", + "t-.t-", ] assert sorted(str(k) for k in a.op_members.keys()) == sorted( [*triv_keys, *keys3] @@ -63,36 +66,8 @@ def test_split_ad_to_evol_map(self): a.op_members[member.MemberName("V.V")].value, ome[(200, 200)].value, ) - # # if alpha is zero, nothing non-trivial should happen - b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, []) - assert sorted(str(k) for k in b.op_members.keys()) == sorted( - [*triv_keys, *keys3] - ) - # assert_almost_equal( - # b.op_members[member.MemberName("V.V")].value, - # np.eye(self.shape[0]), - # ) - # nf=3 + IC - self.update_intrinsic_OME(ome) - c = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [4]) - assert sorted(str(k) for k in c.op_members.keys()) == sorted( - [*triv_keys, *keys3, "S.c+", "g.c+", "c+.c+", "c-.c-"] - ) - assert_almost_equal( - c.op_members[member.MemberName("V.V")].value, - b.op_members[member.MemberName("V.V")].value, - ) - # nf=3 + IB - d = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [5]) - assert sorted(str(k) for k in d.op_members.keys()) == sorted( - [*triv_keys, *keys3, "b+.b+", "b-.b-"] - ) - assert_almost_equal( - d.op_members[member.MemberName("b+.b+")].value, - np.eye(self.shape[0]), - ) # nf=4 + IB - d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1, [5]) + d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1) assert sorted(str(k) for k in d.op_members.keys()) == sorted( [ *triv_keys, @@ -106,6 +81,8 @@ def test_split_ad_to_evol_map(self): "b+.b+", # "b-.V", "b-.b-", + "t+.t+", + "t-.t-", ] ) assert_almost_equal( @@ -127,7 +104,7 @@ def test_split_ad_to_evol_map(self): def test_split_ad_to_evol_map_qed(self): ome = self.mkOME() - a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [], qed=True) + a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, qed=True) triv_keys = [ "ph.ph", "S.S", @@ -144,7 +121,15 @@ def test_split_ad_to_evol_map_qed(self): keys3 = [ "c+.S", "c+.g", + "S.c+", + "g.c+", + "c+.c+", + "c-.c-", # "c-.V", + "b+.b+", + "b-.b-", + "t+.t+", + "t-.t-", ] assert sorted(str(k) for k in a.op_members.keys()) == sorted( [*triv_keys, *keys3] @@ -153,36 +138,8 @@ def test_split_ad_to_evol_map_qed(self): a.op_members[member.MemberName("V.V")].value, ome[(200, 200)].value, ) - # # if alpha is zero, nothing non-trivial should happen - b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [], qed=True) - assert sorted(str(k) for k in b.op_members.keys()) == sorted( - [*triv_keys, *keys3] - ) - # assert_almost_equal( - # b.op_members[member.MemberName("V.V")].value, - # np.eye(self.shape[0]), - # ) - # nf=3 + IC - self.update_intrinsic_OME(ome) - c = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [4], qed=True) - assert sorted(str(k) for k in c.op_members.keys()) == sorted( - [*triv_keys, *keys3, "S.c+", "g.c+", "c+.c+", "c-.c-"] - ) - assert_almost_equal( - c.op_members[member.MemberName("V.V")].value, - b.op_members[member.MemberName("V.V")].value, - ) - # nf=3 + IB - d = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [5], qed=True) - assert sorted(str(k) for k in d.op_members.keys()) == sorted( - [*triv_keys, *keys3, "b+.b+", "b-.b-"] - ) - assert_almost_equal( - d.op_members[member.MemberName("b+.b+")].value, - np.eye(self.shape[0]), - ) # nf=4 + IB - d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1, [5], qed=True) + d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1, qed=True) assert sorted(str(k) for k in d.op_members.keys()) == sorted( [ *triv_keys, @@ -196,6 +153,8 @@ def test_split_ad_to_evol_map_qed(self): "b+.b+", # "b-.V", "b-.b-", + "t+.t+", + "t-.t-", ] ) assert_almost_equal( diff --git a/tests/eko/evolution_operator/test_physical.py b/tests/eko/evolution_operator/test_physical.py index b4a9936b6..fb74fb617 100644 --- a/tests/eko/evolution_operator/test_physical.py +++ b/tests/eko/evolution_operator/test_physical.py @@ -221,27 +221,21 @@ def mk_op_members(shape=(2, 2), qed=False): return om -def get_ad_to_evol_map(nf, intrinsic_range=None, qed=False): +def get_ad_to_evol_map(nf, qed=False): oms = mk_op_members(qed=qed) - m = PhysicalOperator.ad_to_evol_map(oms, nf, 1, intrinsic_range, qed) + m = PhysicalOperator.ad_to_evol_map(oms, nf, 1, qed) return sorted(map(str, m.op_members.keys())) def test_ad_to_evol_map(): triv_ops = ("S.S", "S.g", "g.S", "g.g", "V.V", "V3.V3", "T3.T3", "V8.V8", "T8.T8") # nf=3 - assert sorted(triv_ops) == get_ad_to_evol_map(3) - # nf=3 + IC - assert sorted([*triv_ops, "c+.c+", "c-.c-"]) == get_ad_to_evol_map(3, [4]) - # nf=3 + IC + IB assert sorted( - [*triv_ops, "c+.c+", "c-.c-", "b+.b+", "b-.b-"] - ) == get_ad_to_evol_map(3, [4, 5]) - # nf=4 + IC(non-existant) + IB - ks = sorted([*triv_ops, "V15.V15", "T15.T15", "b+.b+", "b-.b-"]) - assert ks == get_ad_to_evol_map(4, [4, 5]) - # nf=4 + IB - assert ks == get_ad_to_evol_map(4, [5]) + [*triv_ops, "c+.c+", "c-.c-", "b+.b+", "b-.b-", "t+.t+", "t-.t-"] + ) == get_ad_to_evol_map(3) + # nf=4 + ks = sorted([*triv_ops, "V15.V15", "T15.T15", "b+.b+", "b-.b-", "t+.t+", "t-.t-"]) + assert ks == get_ad_to_evol_map(4) # nf=6 assert sorted( [*triv_ops, "T15.T15", "V15.V15", "T24.T24", "V24.V24", "T35.T35", "V35.V35"] @@ -274,18 +268,12 @@ def test_ad_to_evol_map_qed(): "Td3.Td3", ) # nf=3 - assert sorted(triv_ops) == get_ad_to_evol_map(3, qed=True) - # nf=3 + IC - assert sorted([*triv_ops, "c+.c+", "c-.c-"]) == get_ad_to_evol_map(3, [4], qed=True) - # nf=3 + IC + IB assert sorted( - [*triv_ops, "c+.c+", "c-.c-", "b+.b+", "b-.b-"] - ) == get_ad_to_evol_map(3, [4, 5], qed=True) - # nf=4 + IC(non-existant) + IB - ks = sorted([*triv_ops, "Vu3.Vu3", "Tu3.Tu3", "b+.b+", "b-.b-"]) - assert ks == get_ad_to_evol_map(4, [4, 5], qed=True) - # nf=4 + IB - assert ks == get_ad_to_evol_map(4, [5], qed=True) + [*triv_ops, "c+.c+", "c-.c-", "b+.b+", "b-.b-", "t+.t+", "t-.t-"] + ) == get_ad_to_evol_map(3, True) + # nf=4 + ks = sorted([*triv_ops, "Vu3.Vu3", "Tu3.Tu3", "b+.b+", "b-.b-", "t+.t+", "t-.t-"]) + assert ks == get_ad_to_evol_map(4, True) # nf=6 assert sorted( [*triv_ops, "Tu3.Tu3", "Vu3.Vu3", "Td8.Td8", "Vd8.Vd8", "Tu8.Tu8", "Vu8.Vu8"] diff --git a/tests/eko/io/test_runcards.py b/tests/eko/io/test_runcards.py index 6af37edcb..0a391047b 100644 --- a/tests/eko/io/test_runcards.py +++ b/tests/eko/io/test_runcards.py @@ -100,7 +100,6 @@ def test_runcards_quarkmass(): tc["IC"] = 1 oc = copy.deepcopy(operator_card) nt, no = rc.update(tc, oc) - assert nt.heavy.intrinsic_flavors == [4] for _, scale in nt.heavy.masses: assert np.isnan(scale) m2s = rc.masses(nt, no.configs.evolution_method) diff --git a/tests/eko/test_quantities.py b/tests/eko/test_quantities.py index 8fd6c298f..f9dcb9203 100644 --- a/tests/eko/test_quantities.py +++ b/tests/eko/test_quantities.py @@ -34,7 +34,6 @@ def test_HeavyQuarks(): def test_HeavyInfo(): i = hq.HeavyInfo( num_flavs_init=4, - intrinsic_flavors=[4, 5], masses=hq.HeavyQuarkMasses( [ hq.QuarkMassRef([2.0, nan]), From 7279e86333bf09043c8c161144561f6074b3c230 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 1 Aug 2023 11:40:59 +0200 Subject: [PATCH 04/13] Compute always intrinsic in ome --- .../operator_matrix_element.py | 24 +++++++++---------- tests/eko/evolution_operator/test_ome.py | 1 + 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/eko/evolution_operator/operator_matrix_element.py b/src/eko/evolution_operator/operator_matrix_element.py index b45e06c81..fc0518aa1 100644 --- a/src/eko/evolution_operator/operator_matrix_element.py +++ b/src/eko/evolution_operator/operator_matrix_element.py @@ -234,11 +234,10 @@ def labels(self): logger.warning("%s: skipping non-singlet sector", self.log_label) else: labels.append((200, 200)) - if self.backward_method is not None: - # intrinsic labels, which are not zero at NLO - labels.append((br.matching_hminus_pid, br.matching_hminus_pid)) - # These contributions are always 0 for the moment - # labels.extend([(200, br.matching_hminus_pid), (br.matching_hminus_pid, 200)]) + # intrinsic labels, which are not zero at NLO + labels.append((br.matching_hminus_pid, br.matching_hminus_pid)) + # These contributions are always 0 for the moment + # labels.extend([(200, br.matching_hminus_pid), (br.matching_hminus_pid, 200)]) # same for singlet if self.config["debug_skip_singlet"]: logger.warning("%s: skipping singlet sector", self.log_label) @@ -250,14 +249,13 @@ def labels(self): (br.matching_hplus_pid, 100), ] ) - if self.backward_method is not None: - labels.extend( - [ - (21, br.matching_hplus_pid), - (100, br.matching_hplus_pid), - (br.matching_hplus_pid, br.matching_hplus_pid), - ] - ) + labels.extend( + [ + (21, br.matching_hplus_pid), + (100, br.matching_hplus_pid), + (br.matching_hplus_pid, br.matching_hplus_pid), + ] + ) return labels def quad_ker(self, label, logx, areas): diff --git a/tests/eko/evolution_operator/test_ome.py b/tests/eko/evolution_operator/test_ome.py index b844290c9..c90563269 100644 --- a/tests/eko/evolution_operator/test_ome.py +++ b/tests/eko/evolution_operator/test_ome.py @@ -338,6 +338,7 @@ def test_labels(self, theory_ffns, operator_card, tmp_path: pathlib.Path): path = tmp_path / "eko.tar" for skip_singlet in [True, False]: for skip_ns in [True, False]: + operator_card.configs.inversion_method = "exact" operator_card.debug.skip_singlet = skip_singlet operator_card.debug.skip_non_singlet = skip_ns path.unlink(missing_ok=True) From da9fbc42da4ee62823f80031efba037a0e9e8642 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 1 Aug 2023 12:08:15 +0200 Subject: [PATCH 05/13] Move eko.io.runcards.update --- src/eko/io/runcards.py | 22 -------- src/eko/runner/legacy.py | 4 +- src/ekomark/data/__init__.py | 28 ++++++++++ tests/eko/io/test_runcards.py | 88 ------------------------------ tests/ekomark/data/__init__.py | 0 tests/ekomark/data/test_init.py | 96 +++++++++++++++++++++++++++++++++ 6 files changed, 127 insertions(+), 111 deletions(-) create mode 100644 tests/ekomark/data/__init__.py create mode 100644 tests/ekomark/data/test_init.py diff --git a/src/eko/io/runcards.py b/src/eko/io/runcards.py index d5c8c36e8..8a6261bf2 100644 --- a/src/eko/io/runcards.py +++ b/src/eko/io/runcards.py @@ -255,28 +255,6 @@ def new_operator(self): return OperatorCard.from_dict(new) -def update(theory: Union[RawCard, TheoryCard], operator: Union[RawCard, OperatorCard]): - """Update legacy runcards. - - This function is mainly defined for compatibility with the old interface. - Prefer direct usage of :class:`Legacy` in new code. - - Consecutive applications of this function yield identical results:: - - cards = update(theory, operator) - assert update(*cards) == cards - - """ - if isinstance(theory, TheoryCard) or isinstance(operator, OperatorCard): - # if one is not a dict, both have to be new cards - assert isinstance(theory, TheoryCard) - assert isinstance(operator, OperatorCard) - return theory, operator - - cards = Legacy(theory, operator) - return cards.new_theory, cards.new_operator - - def default_atlas(masses: list, matching_ratios: list): r"""Create default landscape. diff --git a/src/eko/runner/legacy.py b/src/eko/runner/legacy.py index 76d6f1a8a..5de53ceb7 100644 --- a/src/eko/runner/legacy.py +++ b/src/eko/runner/legacy.py @@ -3,6 +3,8 @@ import os from typing import Union +from ekomark.data import update_runcards + from ..evolution_operator.grid import OperatorGrid from ..io import EKO, Operator, runcards from ..io.types import RawCard @@ -43,7 +45,7 @@ def __init__( path where to store the computed operator """ - new_theory, new_operator = runcards.update(theory_card, operators_card) + new_theory, new_operator = update_runcards(theory_card, operators_card) # Store inputs self.path = path diff --git a/src/ekomark/data/__init__.py b/src/ekomark/data/__init__.py index 797409204..b6242845b 100644 --- a/src/ekomark/data/__init__.py +++ b/src/ekomark/data/__init__.py @@ -1 +1,29 @@ """EKO database configuration.""" +from typing import Union + +from eko.io.runcards import Legacy, OperatorCard, TheoryCard +from eko.io.types import RawCard + + +def update_runcards( + theory: Union[RawCard, TheoryCard], operator: Union[RawCard, OperatorCard] +): + """Update legacy runcards. + + This function is mainly defined for compatibility with the old interface. + Prefer direct usage of :class:`Legacy` in new code. + + Consecutive applications of this function yield identical results:: + + cards = update(theory, operator) + assert update(*cards) == cards + + """ + if isinstance(theory, TheoryCard) or isinstance(operator, OperatorCard): + # if one is not a dict, both have to be new cards + assert isinstance(theory, TheoryCard) + assert isinstance(operator, OperatorCard) + return theory, operator + + cards = Legacy(theory, operator) + return cards.new_theory, cards.new_operator diff --git a/tests/eko/io/test_runcards.py b/tests/eko/io/test_runcards.py index 0a391047b..465a1abf6 100644 --- a/tests/eko/io/test_runcards.py +++ b/tests/eko/io/test_runcards.py @@ -1,14 +1,11 @@ -import copy from io import StringIO import numpy as np import pytest import yaml -from banana.data.theories import default_card as theory_card from eko.io import runcards as rc from eko.io.bases import Bases -from ekomark.data.operators import default_card as operator_card def test_flavored_mu2grid(): @@ -31,48 +28,6 @@ def test_flavored_mu2grid(): assert t == tuple(l) -def test_runcards_opcard(): - # check conversion - tc = copy.deepcopy(theory_card) - oc = copy.deepcopy(operator_card) - tc["Q0"] = 2.0 - # mugrid - oc["mugrid"] = [2.0, 10.0] - _nt, no = rc.update(tc, oc) - assert isinstance(no, rc.OperatorCard) - assert len(no.evolgrid) == len(oc["mugrid"]) - assert len(no.mu2grid) == len(no.evolgrid) - assert no.evolgrid[0][-1] == 4 - assert no.evolgrid[1][-1] == 5 - np.testing.assert_allclose(no.mu0, tc["Q0"]) - np.testing.assert_allclose(no.mu20, tc["Q0"] ** 2.0) - assert len(no.pids) == 14 - check_dumpable(no) - del oc["mugrid"] - # or mu2grid - oc["mu2grid"] = [9.0, 30.0, 32.0] - _nt, no = rc.update(tc, oc) - assert isinstance(no, rc.OperatorCard) - assert len(no.evolgrid) == len(oc["mu2grid"]) - assert len(no.mu2grid) == len(no.evolgrid) - assert no.evolgrid[0][-1] == 4 - assert no.evolgrid[1][-1] == 5 - assert no.evolgrid[2][-1] == 5 - check_dumpable(no) - del oc["mu2grid"] - # or Q2grid - oc["Q2grid"] = [15.0, 130.0, 140.0, 1e5] - _nt, no = rc.update(tc, oc) - assert isinstance(no, rc.OperatorCard) - assert len(no.evolgrid) == len(oc["Q2grid"]) - assert len(no.mu2grid) == len(no.evolgrid) - assert no.evolgrid[0][-1] == 4 - assert no.evolgrid[1][-1] == 5 - assert no.evolgrid[2][-1] == 5 - assert no.evolgrid[3][-1] == 6 - check_dumpable(no) - - def check_dumpable(no): """Check we can write and read to yaml.""" so = StringIO() @@ -81,49 +36,6 @@ def check_dumpable(no): noo = yaml.safe_load(so) -def test_runcards_ekomark(): - # check conversion - tc = copy.deepcopy(theory_card) - oc = copy.deepcopy(operator_card) - nt, no = rc.update(tc, oc) - assert isinstance(nt, rc.TheoryCard) - assert isinstance(no, rc.OperatorCard) - # check is idempotent - nnt, nno = rc.update(nt, no) - assert nnt == nt - assert nno == no - - -def test_runcards_quarkmass(): - tc = copy.deepcopy(theory_card) - tc["nfref"] = 5 - tc["IC"] = 1 - oc = copy.deepcopy(operator_card) - nt, no = rc.update(tc, oc) - for _, scale in nt.heavy.masses: - assert np.isnan(scale) - m2s = rc.masses(nt, no.configs.evolution_method) - raw = rc.Legacy.heavies("m%s", tc) - raw2 = np.power(raw, 2.0) - np.testing.assert_allclose(m2s, raw2) - tc["HQ"] = "MSBAR" - tc["Qmc"] = raw[0] * 1.1 - tc["Qmb"] = raw[1] * 1.1 - tc["Qmt"] = raw[2] * 0.9 - nt, no = rc.update(tc, oc) - for _, scale in nt.heavy.masses: - assert not np.isnan(scale) - m2s = rc.masses(nt, no.configs.evolution_method) - for m1, m2 in zip(m2s, raw2): - assert not np.isclose(m1, m2) - tc["HQ"] = "Blub" - with pytest.raises(ValueError, match="mass scheme"): - _nt, _no = rc.update(tc, oc) - nt.heavy.masses_scheme = "Bla" - with pytest.raises(ValueError, match="mass scheme"): - _ms = rc.masses(nt, no.configs.evolution_method) - - def test_legacy_fallback(): assert rc.Legacy.fallback(1, 2, 3) == 1 assert rc.Legacy.fallback(None, 2, 3) == 2 diff --git a/tests/ekomark/data/__init__.py b/tests/ekomark/data/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/ekomark/data/test_init.py b/tests/ekomark/data/test_init.py new file mode 100644 index 000000000..7d39dca43 --- /dev/null +++ b/tests/ekomark/data/test_init.py @@ -0,0 +1,96 @@ +import copy + +import numpy as np +import pytest +from banana.data.theories import default_card as theory_card + +from eko.io import runcards as rc +from ekomark.data import update_runcards +from ekomark.data.operators import default_card as operator_card + +from ...eko.io.test_runcards import check_dumpable + + +def test_runcards_opcard(): + # check conversion + tc = copy.deepcopy(theory_card) + oc = copy.deepcopy(operator_card) + tc["Q0"] = 2.0 + # mugrid + oc["mugrid"] = [2.0, 10.0] + _nt, no = update_runcards(tc, oc) + assert isinstance(no, rc.OperatorCard) + assert len(no.evolgrid) == len(oc["mugrid"]) + assert len(no.mu2grid) == len(no.evolgrid) + assert no.evolgrid[0][-1] == 4 + assert no.evolgrid[1][-1] == 5 + np.testing.assert_allclose(no.mu0, tc["Q0"]) + np.testing.assert_allclose(no.mu20, tc["Q0"] ** 2.0) + assert len(no.pids) == 14 + check_dumpable(no) + del oc["mugrid"] + # or mu2grid + oc["mu2grid"] = [9.0, 30.0, 32.0] + _nt, no = update_runcards(tc, oc) + assert isinstance(no, rc.OperatorCard) + assert len(no.evolgrid) == len(oc["mu2grid"]) + assert len(no.mu2grid) == len(no.evolgrid) + assert no.evolgrid[0][-1] == 4 + assert no.evolgrid[1][-1] == 5 + assert no.evolgrid[2][-1] == 5 + check_dumpable(no) + del oc["mu2grid"] + # or Q2grid + oc["Q2grid"] = [15.0, 130.0, 140.0, 1e5] + _nt, no = update_runcards(tc, oc) + assert isinstance(no, rc.OperatorCard) + assert len(no.evolgrid) == len(oc["Q2grid"]) + assert len(no.mu2grid) == len(no.evolgrid) + assert no.evolgrid[0][-1] == 4 + assert no.evolgrid[1][-1] == 5 + assert no.evolgrid[2][-1] == 5 + assert no.evolgrid[3][-1] == 6 + check_dumpable(no) + + +def test_runcards_ekomark(): + # check conversion + tc = copy.deepcopy(theory_card) + oc = copy.deepcopy(operator_card) + nt, no = update_runcards(tc, oc) + assert isinstance(nt, rc.TheoryCard) + assert isinstance(no, rc.OperatorCard) + # check is idempotent + nnt, nno = update_runcards(nt, no) + assert nnt == nt + assert nno == no + + +def test_runcards_quarkmass(): + tc = copy.deepcopy(theory_card) + tc["nfref"] = 5 + tc["IC"] = 1 + oc = copy.deepcopy(operator_card) + nt, no = update_runcards(tc, oc) + for _, scale in nt.heavy.masses: + assert np.isnan(scale) + m2s = rc.masses(nt, no.configs.evolution_method) + raw = rc.Legacy.heavies("m%s", tc) + raw2 = np.power(raw, 2.0) + np.testing.assert_allclose(m2s, raw2) + tc["HQ"] = "MSBAR" + tc["Qmc"] = raw[0] * 1.1 + tc["Qmb"] = raw[1] * 1.1 + tc["Qmt"] = raw[2] * 0.9 + nt, no = update_runcards(tc, oc) + for _, scale in nt.heavy.masses: + assert not np.isnan(scale) + m2s = rc.masses(nt, no.configs.evolution_method) + for m1, m2 in zip(m2s, raw2): + assert not np.isclose(m1, m2) + tc["HQ"] = "Blub" + with pytest.raises(ValueError, match="mass scheme"): + _nt, _no = update_runcards(tc, oc) + nt.heavy.masses_scheme = "Bla" + with pytest.raises(ValueError, match="mass scheme"): + _ms = rc.masses(nt, no.configs.evolution_method) From ac406cd3c0f956654a5f679e95fec0073bb55c00 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 1 Aug 2023 12:40:16 +0200 Subject: [PATCH 06/13] Use EPoint in couplings --- benchmarks/eko/benchmark_evol_to_unity.py | 3 +-- benchmarks/eko/benchmark_msbar_evolution.py | 3 +-- benchmarks/eko/benchmark_strong_coupling.py | 8 +++----- benchmarks/ekobox/benchmark_evol_pdf.py | 4 ++-- doc/source/overview/tutorials/alpha_s.ipynb | 2 +- doc/source/overview/tutorials/pdf.ipynb | 3 +-- extras/n3lo_bench/splitting_function_utils.py | 3 +-- src/eko/couplings.py | 6 +++--- src/eko/io/runcards.py | 3 +-- src/eko/msbar_masses.py | 4 ++-- src/eko/quantities/couplings.py | 10 ++-------- src/ekobox/cards.py | 3 +-- src/ekobox/info_file.py | 2 +- tests/eko/kernels/test_kernels_QEDns.py | 3 +-- tests/eko/quantities/test_couplings.py | 4 ++-- tests/eko/runner/test_legacy.py | 2 +- tests/eko/test_couplings.py | 17 ++++++----------- tests/eko/test_msbar_masses.py | 3 +-- 18 files changed, 31 insertions(+), 52 deletions(-) diff --git a/benchmarks/eko/benchmark_evol_to_unity.py b/benchmarks/eko/benchmark_evol_to_unity.py index d269db202..5e864b3b5 100644 --- a/benchmarks/eko/benchmark_evol_to_unity.py +++ b/benchmarks/eko/benchmark_evol_to_unity.py @@ -19,8 +19,7 @@ def update_cards(theory: TheoryCard, operator: OperatorCard): theory.couplings = CouplingsInfo( alphas=0.35, alphaem=0.007496, - scale=float(np.sqrt(2)), - num_flavs_ref=None, + scale=(float(np.sqrt(2)), 4), ) theory.heavy.num_flavs_init = 4 theory.heavy.masses.c.value = 1.0 diff --git a/benchmarks/eko/benchmark_msbar_evolution.py b/benchmarks/eko/benchmark_msbar_evolution.py index 6286e15ad..86e0e361e 100644 --- a/benchmarks/eko/benchmark_msbar_evolution.py +++ b/benchmarks/eko/benchmark_msbar_evolution.py @@ -20,9 +20,8 @@ def update_theory(theory: TheoryCard): theory.order = (3, 0) - theory.couplings.scale = 91 + theory.couplings.ref = (91, 5) theory.couplings.alphaem = 0.007496 - theory.couplings.num_flavs_ref = 5 theory.heavy.masses_scheme = QuarkMassScheme.MSBAR theory.heavy.masses.c = QuarkMassRef([1.5, 18]) theory.heavy.masses.b = QuarkMassRef([4.1, 20]) diff --git a/benchmarks/eko/benchmark_strong_coupling.py b/benchmarks/eko/benchmark_strong_coupling.py index 80f65cf72..3f94bd036 100644 --- a/benchmarks/eko/benchmark_strong_coupling.py +++ b/benchmarks/eko/benchmark_strong_coupling.py @@ -39,8 +39,7 @@ def ref_couplings(ref_values, ref_scale: float, ref_nf: FlavorsNumber) -> Coupli return CouplingsInfo( alphas=ref_values[0], alphaem=ref_values[1], - scale=ref_scale, - num_flavs_ref=ref_nf, + ref=(ref_scale, ref_nf), ) @@ -841,9 +840,8 @@ def benchmark_APFEL_fact_to_ren_lha_settings(self, theory_card: TheoryCard): theory = theory_card theory.order = (3, 0) theory.couplings.alphas = 0.35 - theory.couplings.scale = float(np.sqrt(2)) theory.couplings.alphaem = 0.007496 - theory.couplings.num_flavs_ref = 4 + theory.couplings.ref = (float(np.sqrt(2)), 4) theory.heavy.num_flavs_init = 3 theory.xif = np.sqrt(1.0 / 2.0) theory.heavy.masses.c.value = np.sqrt(2.0) @@ -886,7 +884,7 @@ def benchmark_APFEL_fact_to_ren_lha_settings(self, theory_card: TheoryCard): apfel.SetTheory("QCD") apfel.SetPerturbativeOrder(theory.order[0] - 1) apfel.SetAlphaEvolution("exact") - apfel.SetAlphaQCDRef(theory.couplings.alphas, theory.couplings.scale) + apfel.SetAlphaQCDRef(theory.couplings.alphas, theory.couplings.ref[0]) apfel.SetVFNS() apfel.SetPoleMasses(qmasses.c.value, qmasses.b.value, qmasses.t.value) apfel.SetMassMatchingScales( diff --git a/benchmarks/ekobox/benchmark_evol_pdf.py b/benchmarks/ekobox/benchmark_evol_pdf.py index c85083604..5d198b3ff 100644 --- a/benchmarks/ekobox/benchmark_evol_pdf.py +++ b/benchmarks/ekobox/benchmark_evol_pdf.py @@ -21,7 +21,7 @@ def benchmark_evolve_single_member( theory = theory_card theory.order = (1, 0) theory.couplings.alphas = 0.118000 - theory.couplings.scale = 91.1876 + theory.couplings.ref = (91.1876, 5) theory.couplings.alphaem = 0.007496 theory.heavy.masses.c.value = 1.3 theory.heavy.masses.b.value = 4.75 @@ -47,7 +47,7 @@ def benchmark_evolve_single_member( ev_pdf = lhapdf.mkPDF("EvPDF", 0) assert info["XMin"] == op.xgrid.raw[0] assert info["SetDesc"] == "MyEvolvedPDF" - assert info["MZ"] == theory.couplings.scale + assert info["MZ"] == theory.couplings.ref[0] assert info["Debug"] == "Debug" xgrid = op.xgrid.raw for idx, mu2 in enumerate(op.mu2grid): diff --git a/doc/source/overview/tutorials/alpha_s.ipynb b/doc/source/overview/tutorials/alpha_s.ipynb index d8954900d..327eb39c9 100644 --- a/doc/source/overview/tutorials/alpha_s.ipynb +++ b/doc/source/overview/tutorials/alpha_s.ipynb @@ -40,7 +40,7 @@ "\n", "# set the (alpha_s, alpha_em) reference values\n", "couplings_ref = CouplingsInfo(\n", - " alphas=0.118, alphaem=0.007496252, scale=91.0, num_flavs_ref=5\n", + " alphas=0.118, alphaem=0.007496252, ref=(91.0, 5)\n", ")\n", "\n", "# set heavy quark masses and their threshold ratios\n", diff --git a/doc/source/overview/tutorials/pdf.ipynb b/doc/source/overview/tutorials/pdf.ipynb index 0719d9d63..2cc2f688b 100644 --- a/doc/source/overview/tutorials/pdf.ipynb +++ b/doc/source/overview/tutorials/pdf.ipynb @@ -369,8 +369,7 @@ "th_card.orders = (1,0) # QCD LO\n", "th_card.heavy.masses = HeavyQuarks([QuarkMassRef([1.3,nan]), QuarkMassRef([4.75,nan]), QuarkMassRef([172.,nan])]) # quark mass\n", "th_card.couplings.alphas = 0.130000 # reference value of alpha_s\n", - "th_card.couplings.scale = 91.1876 # the reference scale at which alpha_s is provided\n", - "th_card.couplings.num_flavs_ref = 5 # the number of flavors active at the alpha_s reference scale\n", + "th_card.couplings.ref = (91.1876,5) # the reference scale together with the number of flavors at which alpha_s is provided\n", "th_card.heavy.num_flavs_init = 3 # the number of flavors active at the reference scale" ] }, diff --git a/extras/n3lo_bench/splitting_function_utils.py b/extras/n3lo_bench/splitting_function_utils.py index db06f5c72..7368e24b1 100644 --- a/extras/n3lo_bench/splitting_function_utils.py +++ b/extras/n3lo_bench/splitting_function_utils.py @@ -103,8 +103,7 @@ def compute_a_s(q2=None, xif2=1.0, nf=None, order=(4, 0)): ref = CouplingsInfo( alphas=0.1181, alphaem=0.007496, - scale=91.00, - num_flavs_ref=5, + ref=(91.00, 5), ) sc = Couplings( couplings=ref, diff --git a/src/eko/couplings.py b/src/eko/couplings.py index 1b3620032..3bf21b11d 100644 --- a/src/eko/couplings.py +++ b/src/eko/couplings.py @@ -435,7 +435,7 @@ def assert_positive(name, var): assert_positive("alpha_s_ref", couplings.alphas) assert_positive("alpha_em_ref", couplings.alphaem) - assert_positive("scale_ref", couplings.scale) + assert_positive("scale_ref", couplings.ref[0]) if order[0] not in [1, 2, 3, 4]: raise NotImplementedError("a_s beyond N3LO is not implemented") if order[1] not in [0, 1, 2]: @@ -445,7 +445,7 @@ def assert_positive(name, var): raise ValueError(f"Unknown method {method.value}") self.method = method.value - nf_ref = couplings.num_flavs_ref + nf_ref = couplings.ref[1] scheme_name = hqm_scheme.name self.alphaem_running = couplings.em_running self.decoupled_running = False @@ -454,7 +454,7 @@ def assert_positive(name, var): self.a_ref = np.array(couplings.values) / 4.0 / np.pi # convert to a_s and a_em matching_scales = (np.array(masses) * np.array(thresholds_ratios)).tolist() self.thresholds_ratios = list(thresholds_ratios) - self.atlas = matchings.Atlas(matching_scales, (couplings.scale**2, nf_ref)) + self.atlas = matchings.Atlas(matching_scales, (couplings.ref[0] ** 2, nf_ref)) self.hqm_scheme = scheme_name logger.info( "Strong Coupling: a_s(µ_R^2=%f)%s=%f=%f/(4π)", diff --git a/src/eko/io/runcards.py b/src/eko/io/runcards.py index 8a6261bf2..3d8413d1b 100644 --- a/src/eko/io/runcards.py +++ b/src/eko/io/runcards.py @@ -183,8 +183,7 @@ def new_theory(self): alphas=old["alphas"], alphaem=alphaem, em_running=em_running, - scale=old["Qref"], - num_flavs_ref=old["nfref"], + ref=(old["Qref"], old["nfref"]), ) new["heavy"] = { "num_flavs_init": nf_default(old["Q0"] ** 2.0, default_atlas(ms, ks)) diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index f13582110..5edac4c46 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -377,8 +377,8 @@ def compute( """ # TODO: sketch in the docs how the MSbar computation works with a figure. - mu2_ref = couplings.scale**2 - nf_ref: FlavorsNumber = couplings.num_flavs_ref + mu2_ref = couplings.ref[0] ** 2 + nf_ref: FlavorsNumber = couplings.ref[1] masses = np.concatenate((np.zeros(nf_ref - 3), np.full(6 - nf_ref, np.inf))) def sc(thr_masses): diff --git a/src/eko/quantities/couplings.py b/src/eko/quantities/couplings.py index f80420a9c..3a5ab96ef 100644 --- a/src/eko/quantities/couplings.py +++ b/src/eko/quantities/couplings.py @@ -3,6 +3,7 @@ import enum from ..io.dictlike import DictLike +from ..io.types import EvolutionPoint as EPoint from ..io.types import FlavorsNumber, LinearScale, ReferenceRunning, Scalar Coupling = Scalar @@ -19,14 +20,7 @@ class CouplingsInfo(DictLike): alphas: Coupling alphaem: Coupling - scale: LinearScale - num_flavs_ref: FlavorsNumber - r"""Number of active flavors at strong coupling reference scale. - - I.e. :math:`n_{f,\text{ref}}(\mu^2_{\text{ref}})`, formerly called - ``nfref``. - - """ + ref: EPoint em_running: bool = False @property diff --git a/src/ekobox/cards.py b/src/ekobox/cards.py index 482e65faf..3d4e346aa 100644 --- a/src/ekobox/cards.py +++ b/src/ekobox/cards.py @@ -13,8 +13,7 @@ couplings=dict( alphas=0.118, alphaem=0.007496252, - scale=91.2, - num_flavs_ref=5, + ref=(91.2, 5), ), heavy=dict( num_flavs_init=4, diff --git a/src/ekobox/info_file.py b/src/ekobox/info_file.py index 2f99d8377..35f2cfcae 100644 --- a/src/ekobox/info_file.py +++ b/src/ekobox/info_file.py @@ -48,7 +48,7 @@ def build( template_info["OrderQCD"] = theory_card.order[0] - 1 template_info["QMin"] = round(math.sqrt(operators_card.mu2grid[0]), 4) template_info["QMax"] = round(math.sqrt(operators_card.mu2grid[-1]), 4) - template_info["MZ"] = theory_card.couplings.scale + template_info["MZ"] = theory_card.couplings.ref[0] template_info["MUp"] = 0.0 template_info["MDown"] = 0.0 template_info["MStrange"] = 0.0 diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py index b85a3f013..2e710c722 100644 --- a/tests/eko/kernels/test_kernels_QEDns.py +++ b/tests/eko/kernels/test_kernels_QEDns.py @@ -119,8 +119,7 @@ def test_zero_true_gamma(): dict( alphas=alpharef[0], alphaem=alpharef[1], - scale=muref, - num_flavs_ref=5, + ref=(muref, 5), ) ) evmod = CouplingEvolutionMethod.EXACT diff --git a/tests/eko/quantities/test_couplings.py b/tests/eko/quantities/test_couplings.py index 83bef709c..e4d194d5b 100644 --- a/tests/eko/quantities/test_couplings.py +++ b/tests/eko/quantities/test_couplings.py @@ -3,7 +3,7 @@ def test_couplings_ref(): scale = 90.0 - d = dict(alphas=0.1, alphaem=0.01, scale=scale, num_flavs_ref=5) + d = dict(alphas=0.1, alphaem=0.01, ref=(scale, 5)) couplings = CouplingsInfo.from_dict(d) - assert couplings.scale == scale + assert couplings.ref[0] == scale assert not couplings.em_running diff --git a/tests/eko/runner/test_legacy.py b/tests/eko/runner/test_legacy.py index 60395acc5..1828285c3 100644 --- a/tests/eko/runner/test_legacy.py +++ b/tests/eko/runner/test_legacy.py @@ -36,7 +36,7 @@ class FakeEM(enum.Enum): eko.runner.legacy.Runner(theory_card, operator_card, path=path) # MSbar scheme theory_card.heavy.masses_scheme = QuarkMassScheme.MSBAR - theory_card.couplings.num_flavs_ref = 5 + theory_card.couplings.ref = (91.0, 5) theory_card.heavy.masses.c.scale = 2 theory_card.heavy.masses.b.scale = 4.5 theory_card.heavy.masses.t.scale = 173.07 diff --git a/tests/eko/test_couplings.py b/tests/eko/test_couplings.py index 10535cbcb..7468051be 100644 --- a/tests/eko/test_couplings.py +++ b/tests/eko/test_couplings.py @@ -53,8 +53,7 @@ def test_init(self): dict( alphas=alpharef[0], alphaem=alpharef[1], - scale=muref, - num_flavs_ref=5, + ref=(muref, 5), ) ) order = (1, 0) @@ -107,7 +106,7 @@ def test_init(self): ) with pytest.raises(ValueError): coup3 = copy.deepcopy(couplings) - coup3.scale = 0 + coup3.ref = (0.0, 5) Couplings( coup3, order, @@ -159,8 +158,7 @@ def test_ref(self): dict( alphas=alpharef[0], alphaem=alpharef[1], - scale=muref, - num_flavs_ref=nfref, + ref=(muref, nfref), ) ) for order_s in [1, 2, 3, 4]: @@ -193,8 +191,7 @@ def test_ref_copy_e841b0dfdee2f31d9ccc1ecee4d9d1a6f6624313(self): dict( alphas=alpharef[0], alphaem=alpharef[1], - scale=muref, - num_flavs_ref=3, # reference nf is needed to force the matching + ref=(muref, 3), # reference nf is needed to force the matching ) ) sc = Couplings( @@ -233,8 +230,7 @@ def test_exact(self): dict( alphas=alpharef[0], alphaem=alpharef[1], - scale=muref, - num_flavs_ref=nfref, + ref=(muref, nfref), em_running=em_running, ) ) @@ -295,8 +291,7 @@ def benchmark_expanded_n3lo(self): couplings = CouplingsInfo( alphas=alpharef[0], alphaem=alpharef[1], - scale=muref, - num_flavs_ref=5, + ref=(muref, 5), ) m2c = 2 m2b = 25 diff --git a/tests/eko/test_msbar_masses.py b/tests/eko/test_msbar_masses.py index b0c102ef8..edfd243ab 100644 --- a/tests/eko/test_msbar_masses.py +++ b/tests/eko/test_msbar_masses.py @@ -15,9 +15,8 @@ def theory_card(theory_card: TheoryCard): th = theory_card th.order = (3, 0) th.couplings.alphas = 0.1180 - th.couplings.scale = 91.0 th.couplings.alphaem = 0.00781 - th.couplings.num_flavs_ref = 5 + th.couplings.ref = (91.0, 5) th.heavy.masses = HeavyQuarkMasses( [QuarkMassRef(val) for val in [(2.0, 2.1), (4.0, 4.1), (175.0, 174.9)]] ) From 45ae33043c33435ae5231e34785d1572cd774fe7 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 1 Aug 2023 15:04:11 +0200 Subject: [PATCH 07/13] Fix couplings ref ep in benchmarks --- benchmarks/eko/benchmark_evol_to_unity.py | 2 +- benchmarks/eko/benchmark_msbar_evolution.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/benchmarks/eko/benchmark_evol_to_unity.py b/benchmarks/eko/benchmark_evol_to_unity.py index 5e864b3b5..561721f6f 100644 --- a/benchmarks/eko/benchmark_evol_to_unity.py +++ b/benchmarks/eko/benchmark_evol_to_unity.py @@ -19,7 +19,7 @@ def update_cards(theory: TheoryCard, operator: OperatorCard): theory.couplings = CouplingsInfo( alphas=0.35, alphaem=0.007496, - scale=(float(np.sqrt(2)), 4), + ref=(float(np.sqrt(2)), 4), ) theory.heavy.num_flavs_init = 4 theory.heavy.masses.c.value = 1.0 diff --git a/benchmarks/eko/benchmark_msbar_evolution.py b/benchmarks/eko/benchmark_msbar_evolution.py index 86e0e361e..c0e849e64 100644 --- a/benchmarks/eko/benchmark_msbar_evolution.py +++ b/benchmarks/eko/benchmark_msbar_evolution.py @@ -148,7 +148,7 @@ def benchmark_APFEL_msbar_evolution( apfel.SetTheory("QCD") apfel.SetPerturbativeOrder(order - 1) apfel.SetAlphaEvolution(method) - apfel.SetAlphaQCDRef(coupl.alphas, coupl.scale) + apfel.SetAlphaQCDRef(coupl.alphas, coupl.ref[0]) apfel.SetVFNS() apfel.SetMSbarMasses( qmasses.c.value, qmasses.b.value, qmasses.t.value @@ -216,7 +216,7 @@ def benchmark_APFEL_msbar_solution( apfel.SetTheory("QCD") apfel.SetPerturbativeOrder(order - 1) apfel.SetAlphaEvolution("exact") - apfel.SetAlphaQCDRef(coupl.alphas, coupl.scale) + apfel.SetAlphaQCDRef(coupl.alphas, coupl.ref[0]) apfel.SetVFNS() apfel.SetMSbarMasses(qmasses.c.value, qmasses.b.value, qmasses.t.value) apfel.SetMassScaleReference( From 0f94cdc66171c769c14dbd77052efdc44ea4f029 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 1 Aug 2023 15:38:24 +0200 Subject: [PATCH 08/13] Use EPoint for PDF boundary --- benchmarks/eko/benchmark_evol_to_unity.py | 3 +-- benchmarks/eko/benchmark_strong_coupling.py | 1 - benchmarks/ekobox/benchmark_evol_pdf.py | 4 ++-- doc/source/overview/tutorials/pdf.ipynb | 5 ++--- src/eko/io/legacy.py | 7 +++---- src/eko/io/runcards.py | 22 ++++++++++----------- src/eko/io/struct.py | 2 +- src/eko/quantities/heavy_quarks.py | 6 ------ src/eko/runner/commons.py | 2 +- src/ekobox/cards.py | 3 +-- src/ekobox/cli/runcards.py | 2 +- src/ekobox/info_file.py | 4 +++- src/ekobox/utils.py | 2 +- tests/eko/evolution_operator/test_grid.py | 2 +- tests/eko/io/test_manipulate.py | 2 +- tests/eko/test_quantities.py | 1 - tests/ekobox/test_cards.py | 8 ++++---- tests/ekobox/test_evol_pdf.py | 2 +- tests/ekobox/test_info_file.py | 2 +- tests/ekobox/test_utils.py | 7 +++---- tests/ekomark/data/test_init.py | 2 +- 21 files changed, 39 insertions(+), 50 deletions(-) diff --git a/benchmarks/eko/benchmark_evol_to_unity.py b/benchmarks/eko/benchmark_evol_to_unity.py index 561721f6f..c466d00d8 100644 --- a/benchmarks/eko/benchmark_evol_to_unity.py +++ b/benchmarks/eko/benchmark_evol_to_unity.py @@ -21,11 +21,10 @@ def update_cards(theory: TheoryCard, operator: OperatorCard): alphaem=0.007496, ref=(float(np.sqrt(2)), 4), ) - theory.heavy.num_flavs_init = 4 theory.heavy.masses.c.value = 1.0 theory.heavy.masses.b.value = 4.75 theory.heavy.masses.t.value = 173.0 - operator.mu0 = float(np.sqrt(2)) + operator.init = (float(np.sqrt(2)), 4) operator.mugrid = [(10, 5)] operator.xgrid = XGrid(np.linspace(1e-1, 1, 30)) operator.configs.interpolation_polynomial_degree = 1 diff --git a/benchmarks/eko/benchmark_strong_coupling.py b/benchmarks/eko/benchmark_strong_coupling.py index 3f94bd036..1d1cd7219 100644 --- a/benchmarks/eko/benchmark_strong_coupling.py +++ b/benchmarks/eko/benchmark_strong_coupling.py @@ -842,7 +842,6 @@ def benchmark_APFEL_fact_to_ren_lha_settings(self, theory_card: TheoryCard): theory.couplings.alphas = 0.35 theory.couplings.alphaem = 0.007496 theory.couplings.ref = (float(np.sqrt(2)), 4) - theory.heavy.num_flavs_init = 3 theory.xif = np.sqrt(1.0 / 2.0) theory.heavy.masses.c.value = np.sqrt(2.0) theory.heavy.masses.b.value = 4.5 diff --git a/benchmarks/ekobox/benchmark_evol_pdf.py b/benchmarks/ekobox/benchmark_evol_pdf.py index 5d198b3ff..a0ddcaac4 100644 --- a/benchmarks/ekobox/benchmark_evol_pdf.py +++ b/benchmarks/ekobox/benchmark_evol_pdf.py @@ -27,7 +27,7 @@ def benchmark_evolve_single_member( theory.heavy.masses.b.value = 4.75 theory.heavy.masses.t.value = 172 op = operator_card - op.mu0 = 5.0 + op.init = (5.0, 4) op.mugrid = mugrid # lhapdf import (maybe i have to dump with a x*), do plots) with lhapdf_path(test_pdf): @@ -74,7 +74,7 @@ def benchmark_evolve_more_members( theory = theory_card theory.order = (1, 0) op = operator_card - op.mu0 = 1.0 + op.init = (1.0, 3) op.mugrid = [(10.0, 5), (100.0, 5)] op.xgrid = XGrid([1e-7, 0.01, 0.1, 0.2, 0.3]) with lhapdf_path(test_pdf): diff --git a/doc/source/overview/tutorials/pdf.ipynb b/doc/source/overview/tutorials/pdf.ipynb index 2cc2f688b..923520f59 100644 --- a/doc/source/overview/tutorials/pdf.ipynb +++ b/doc/source/overview/tutorials/pdf.ipynb @@ -361,7 +361,7 @@ "op_card = example.operator()\n", "op_card.xgrid = eko.interpolation.XGrid(make_grid(30, 30)) # x grid\n", "op_card.mugrid = [(float(q),5) for q in np.geomspace(5., 100, 5)] # Q2 grid\n", - "op_card.mu0 = 1.295000 # starting point for the evolution \n", + "op_card.init = (1.295000,3) # starting point for the evolution \n", "\n", "# setup the theory card - this can be mostly inferred from the PDF's .info file\n", "\n", @@ -369,8 +369,7 @@ "th_card.orders = (1,0) # QCD LO\n", "th_card.heavy.masses = HeavyQuarks([QuarkMassRef([1.3,nan]), QuarkMassRef([4.75,nan]), QuarkMassRef([172.,nan])]) # quark mass\n", "th_card.couplings.alphas = 0.130000 # reference value of alpha_s\n", - "th_card.couplings.ref = (91.1876,5) # the reference scale together with the number of flavors at which alpha_s is provided\n", - "th_card.heavy.num_flavs_init = 3 # the number of flavors active at the reference scale" + "th_card.couplings.ref = (91.1876,5) # the reference scale together with the number of flavors at which alpha_s is provided" ] }, { diff --git a/src/eko/io/legacy.py b/src/eko/io/legacy.py index 8491a6acc..3bee940ac 100644 --- a/src/eko/io/legacy.py +++ b/src/eko/io/legacy.py @@ -99,7 +99,6 @@ class PseudoTheory(DictLike): def from_old(cls, old: RawCard): """Load from old metadata.""" heavy = HeavyInfo( - num_flavs_init=4, masses=HeavyQuarkMasses( [ ReferenceRunning([_MC, np.inf]), @@ -122,7 +121,7 @@ class PseudoOperator(DictLike): """ - mu20: float + init: EPoint evolgrid: List[EPoint] xgrid: XGrid configs: dict @@ -130,7 +129,7 @@ class PseudoOperator(DictLike): @classmethod def from_old(cls, old: RawCard): """Load from old metadata.""" - mu20 = float(old["q2_ref"]) + mu0 = float(np.sqrt(float(old["q2_ref"]))) mu2list = old.get("Q2grid") if mu2list is None: mu2list = old["mu2grid"] @@ -146,7 +145,7 @@ def from_old(cls, old: RawCard): interpolation_is_log=old.get("interpolation_is_log"), ) - return cls(mu20=mu20, evolgrid=evolgrid, xgrid=xgrid, configs=configs) + return cls(init=(mu0, 4), evolgrid=evolgrid, xgrid=xgrid, configs=configs) ARRAY_SUFFIX = ".npy.lz4" diff --git a/src/eko/io/runcards.py b/src/eko/io/runcards.py index 3d8413d1b..8ab9beeff 100644 --- a/src/eko/io/runcards.py +++ b/src/eko/io/runcards.py @@ -95,7 +95,7 @@ class Configs(DictLike): class OperatorCard(DictLike): """Operator Card info.""" - mu0: float + init: EPoint """Initial scale.""" mugrid: List[EPoint] xgrid: interpolation.XGrid @@ -115,7 +115,7 @@ class OperatorCard(DictLike): @property def mu20(self): """Squared value of initial scale.""" - return self.mu0**2 + return self.init[0] ** 2 @property def mu2grid(self) -> npt.NDArray: @@ -186,9 +186,6 @@ def new_theory(self): ref=(old["Qref"], old["nfref"]), ) new["heavy"] = { - "num_flavs_init": nf_default(old["Q0"] ** 2.0, default_atlas(ms, ks)) - if old["nf0"] is None - else old["nf0"], "matching_ratios": self.heavies("k%sThr", old), "masses_scheme": old["HQ"], } @@ -212,17 +209,20 @@ def new_operator(self): old_th = self.theory new = {} - new["mu0"] = old_th["Q0"] + ms = self.heavies("m%s", old_th) + ks = self.heavies("k%sThr", old_th) + new["init"] = ( + old_th["Q0"], + nf_default(old_th["Q0"] ** 2.0, default_atlas(ms, ks)) + if old_th["nf0"] is None + else old_th["nf0"], + ) if "mugrid" in old: mugrid = old["mugrid"] else: mu2grid = old["Q2grid"] if "Q2grid" in old else old["mu2grid"] mugrid = np.sqrt(mu2grid).tolist() - new["mugrid"] = flavored_mugrid( - mugrid, - list(self.heavies("m%s", old_th)), - list(self.heavies("k%sThr", old_th)), - ) + new["mugrid"] = flavored_mugrid(mugrid, list(ms), list(ks)) new["configs"] = {} evmod = old_th["ModEv"] diff --git a/src/eko/io/struct.py b/src/eko/io/struct.py index ce4679871..783e82f51 100644 --- a/src/eko/io/struct.py +++ b/src/eko/io/struct.py @@ -552,7 +552,7 @@ def build(self) -> EKO: self.access.open = True metadata = Metadata( _path=self.path, - origin=(self.operator.mu20, self.theory.heavy.num_flavs_init), + origin=(self.operator.init[0] ** 2, self.operator.init[1]), bases=Bases(xgrid=self.operator.xgrid), ) InternalPaths(self.path).bootstrap( diff --git a/src/eko/quantities/heavy_quarks.py b/src/eko/quantities/heavy_quarks.py index 3b081a03d..fe20fc53c 100644 --- a/src/eko/quantities/heavy_quarks.py +++ b/src/eko/quantities/heavy_quarks.py @@ -84,12 +84,6 @@ class HeavyInfo(DictLike): """ - num_flavs_init: FlavorsNumber - r"""Number of active flavors at fitting scale. - - I.e. :math:`n_{f,\text{ref}}(\mu^2_0)`, formerly called ``nf0``. - - """ masses: HeavyQuarkMasses """List of heavy quark masses.""" masses_scheme: QuarkMassScheme diff --git a/src/eko/runner/commons.py b/src/eko/runner/commons.py index c27b4c8a2..05a1d8dd3 100644 --- a/src/eko/runner/commons.py +++ b/src/eko/runner/commons.py @@ -32,7 +32,7 @@ def atlas(theory: TheoryCard, operator: OperatorCard) -> Atlas: # TODO: cache result masses = runcards.masses(theory, operator.configs.evolution_method) matching_scales = np.power(theory.heavy.matching_ratios, 2.0) * np.array(masses) - return Atlas(matching_scales.tolist(), (operator.mu20, theory.heavy.num_flavs_init)) + return Atlas(matching_scales.tolist(), (operator.mu20, operator.init[1])) def couplings(theory: TheoryCard, operator: OperatorCard) -> Couplings: diff --git a/src/ekobox/cards.py b/src/ekobox/cards.py index 3d4e346aa..d56e50da0 100644 --- a/src/ekobox/cards.py +++ b/src/ekobox/cards.py @@ -16,7 +16,6 @@ ref=(91.2, 5), ), heavy=dict( - num_flavs_init=4, masses=[ReferenceRunning([mq, nan]) for mq in (2.0, 4.5, 173.07)], masses_scheme="POLE", matching_ratios=[1.0, 1.0, 1.0], @@ -26,7 +25,7 @@ ) _operator = dict( - mu0=1.65, + init=(1.65, 4), mugrid=[(100.0, 5)], xgrid=np.geomspace(1e-7, 1.0, 50).tolist(), configs=dict( diff --git a/src/ekobox/cli/runcards.py b/src/ekobox/cli/runcards.py index 7fc46f1ef..c5511add5 100644 --- a/src/ekobox/cli/runcards.py +++ b/src/ekobox/cli/runcards.py @@ -35,7 +35,7 @@ def sub_example(destination: pathlib.Path): theory.order = (1, 0) cards.dump(theory.raw, path=destination / "theory.yaml") operator = cards.example.operator() - operator.mu0 = 1.65 + operator.init = (1.65, 4) operator.mugrid = [(np.sqrt(1e5), 5)] cards.dump(operator.raw, path=destination / "operator.yaml") _logger.info(f"Runcards generated to '{destination}'") diff --git a/src/ekobox/info_file.py b/src/ekobox/info_file.py index 35f2cfcae..3d4ca2634 100644 --- a/src/ekobox/info_file.py +++ b/src/ekobox/info_file.py @@ -35,7 +35,9 @@ def build( info file in lhapdf format """ template_info = copy.deepcopy(load.template_info) - template_info["SetDesc"] = "Evolved PDF from " + str(operators_card.mu0) + " GeV" + template_info["SetDesc"] = ( + "Evolved PDF from " + str(operators_card.init[0]) + " GeV" + ) template_info["Authors"] = "" template_info["FlavorScheme"] = "variable" template_info.update(info_update) diff --git a/src/ekobox/utils.py b/src/ekobox/utils.py index 095e605e9..4dc151e96 100644 --- a/src/ekobox/utils.py +++ b/src/ekobox/utils.py @@ -37,7 +37,7 @@ def ekos_product( # another kind of output which includes the theory and operator runcards) ep_match = eko_ini.approx( - (eko_fin.operator_card.mu0**2, eko_fin.theory_card.heavy.num_flavs_init), + (eko_fin.operator_card.init[0] ** 2, eko_fin.operator_card.init[1]), rtol=rtol, atol=atol, ) diff --git a/tests/eko/evolution_operator/test_grid.py b/tests/eko/evolution_operator/test_grid.py index 804388db6..01754578a 100644 --- a/tests/eko/evolution_operator/test_grid.py +++ b/tests/eko/evolution_operator/test_grid.py @@ -72,7 +72,7 @@ def test_mod_expanded(theory_card, theory_ffns, operator_card, tmp_path: pathlib else: theory = theory_card theory.order = (1, 0) - theory.heavy.num_flavs_init = nf0 + operator_card.init = (operator_card.init[0], nf0) path.unlink(missing_ok=True) opgrid = legacy.Runner(theory, operator_card, path=path).op_grid opg = opgrid.compute() diff --git a/tests/eko/io/test_manipulate.py b/tests/eko/io/test_manipulate.py index 3431d93ce..08b38e23d 100644 --- a/tests/eko/io/test_manipulate.py +++ b/tests/eko/io/test_manipulate.py @@ -213,7 +213,7 @@ def _test_to_all_evol( mu2_out = mu_out**2 nfout = 4 epout = (mu2_out, nfout) - eko_factory.operator.mu0 = float(np.sqrt(1.0)) + eko_factory.operator.init = (float(np.sqrt(1.0)), 3) eko_factory.operator.mugrid = [(mu_out, nfout)] eko_factory.operator.xgrid = xgrid eko_factory.operator.configs.interpolation_polynomial_degree = 1 diff --git a/tests/eko/test_quantities.py b/tests/eko/test_quantities.py index f9dcb9203..d3c412a25 100644 --- a/tests/eko/test_quantities.py +++ b/tests/eko/test_quantities.py @@ -33,7 +33,6 @@ def test_HeavyQuarks(): def test_HeavyInfo(): i = hq.HeavyInfo( - num_flavs_init=4, masses=hq.HeavyQuarkMasses( [ hq.QuarkMassRef([2.0, nan]), diff --git a/tests/ekobox/test_cards.py b/tests/ekobox/test_cards.py index 1a35e137e..00ef0e4c3 100644 --- a/tests/ekobox/test_cards.py +++ b/tests/ekobox/test_cards.py @@ -10,14 +10,14 @@ def test_generate_ocard(): mu0 = 1.65 mugrid = [(10.0, 6), (100.0, 5)] op = cards.example.operator() - op.mu0 = mu0 + op.init = (mu0, 4) op.mugrid = mugrid assert pytest.approx(op.mugrid) == mugrid assert pytest.approx(op.mu2grid) == np.array([mu**2 for mu, _ in mugrid]) assert op.configs.interpolation_polynomial_degree == 4 mugrid1 = [100.0, 5] op = cards.example.operator() - op.mu0 = mu0 + op.init = (mu0, 4) op.mugrid = mugrid1 op.configs.interpolation_polynomial_degree = 2 op.configs.interpolation_is_log = False @@ -35,7 +35,7 @@ def test_dump_load_op_card(tmp_path, cd): path1 = tmp_path / "debug_op.yaml" path2 = tmp_path / "debug_op_two.yaml" op = cards.example.operator() - op.mu0 = mu0 + op.init = (mu0, 4) cards.dump(op.raw, path1) cards.dump(op.raw, path2) op_loaded = cards.load(path1) @@ -52,7 +52,7 @@ def test_generate_theory_card(): assert theory.order[0] == 2 rawt = cards.example.raw_theory() assert isinstance(rawt, dict) - assert theory.heavy.num_flavs_init == rawt["heavy"]["num_flavs_init"] + assert theory.heavy.masses.c[0] == rawt["heavy"]["masses"][0][0] def containsnan(obj) -> bool: diff --git a/tests/ekobox/test_evol_pdf.py b/tests/ekobox/test_evol_pdf.py index 1c51dfdad..8e51c6e43 100644 --- a/tests/ekobox/test_evol_pdf.py +++ b/tests/ekobox/test_evol_pdf.py @@ -9,7 +9,7 @@ def init_cards(): op = cards.example.operator() - op.mu0 = 1.65 + op.init = (1.65, 4) op.xgrid = XGrid([0.1, 0.5, 1.0]) op.configs.interpolation_polynomial_degree = 1 theory = cards.example.theory() diff --git a/tests/ekobox/test_info_file.py b/tests/ekobox/test_info_file.py index 64af8eb93..ac50a2253 100644 --- a/tests/ekobox/test_info_file.py +++ b/tests/ekobox/test_info_file.py @@ -10,7 +10,7 @@ def test_build(): theory.order = (2, 0) theory.couplings.alphas = 0.2 op = cards.example.operator() - op.mu0 = 1.0 + op.init = (1.0, 3) op.mugrid = [(10.0, 5), (100.0, 5)] info = info_file.build( theory, op, 4, info_update={"SetDesc": "Prova", "NewArg": 15.3, "MTop": 1.0} diff --git a/tests/ekobox/test_utils.py b/tests/ekobox/test_utils.py index 557db679c..1233d1d91 100644 --- a/tests/ekobox/test_utils.py +++ b/tests/ekobox/test_utils.py @@ -16,10 +16,9 @@ def test_ekos_product(tmp_path): theory = cards.example.theory() theory.order = (1, 0) - theory.heavy.num_flavs_init = 5 op1 = cards.example.operator() - op1.mu0 = mu01 + op1.init = (mu01, 5) op1.mugrid = mugrid1 op1.xgrid = xgrid op1.configs.interpolation_polynomial_degree = 1 @@ -28,13 +27,13 @@ def test_ekos_product(tmp_path): mugrid2 = [(8.0, 5), (10.0, 5), (12.0, 5)] op2 = cards.example.operator() - op2.mu0 = mu0 + op2.init = (mu0, 5) op2.mugrid = mugrid2 op2.xgrid = xgrid op2.configs.interpolation_polynomial_degree = 1 op_err = copy.deepcopy(op2) - op_err.mu0 = mu01 + op_err.init = (mu01, 5) mu2first = (mugrid2[0][0] ** 2, mugrid2[0][1]) diff --git a/tests/ekomark/data/test_init.py b/tests/ekomark/data/test_init.py index 7d39dca43..278a06181 100644 --- a/tests/ekomark/data/test_init.py +++ b/tests/ekomark/data/test_init.py @@ -24,7 +24,7 @@ def test_runcards_opcard(): assert len(no.mu2grid) == len(no.evolgrid) assert no.evolgrid[0][-1] == 4 assert no.evolgrid[1][-1] == 5 - np.testing.assert_allclose(no.mu0, tc["Q0"]) + np.testing.assert_allclose(no.init[0], tc["Q0"]) np.testing.assert_allclose(no.mu20, tc["Q0"] ** 2.0) assert len(no.pids) == 14 check_dumpable(no) From 0d7a459d01b57fcbe34e0875f2e59a1d31373a76 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 1 Aug 2023 16:00:20 +0200 Subject: [PATCH 09/13] Fix benchmark_APFEL_vfns_fact_to_ren --- benchmarks/eko/benchmark_strong_coupling.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/benchmarks/eko/benchmark_strong_coupling.py b/benchmarks/eko/benchmark_strong_coupling.py index 1d1cd7219..d60bd28b0 100644 --- a/benchmarks/eko/benchmark_strong_coupling.py +++ b/benchmarks/eko/benchmark_strong_coupling.py @@ -334,7 +334,7 @@ def benchmark_APFEL_vfns_fact_to_ren(self): ] coupling_ref = np.array([0.118, 0.007496]) scale_ref = 91.0 - nf_ref = 5 + nf_refs = (5, 6) fact_to_ren_lin_list = [0.567, 2.34] masses_list = np.power([2, 2 * 4, 2 * 92], 2) apfel_vals_dict_list = [ @@ -427,8 +427,8 @@ def benchmark_APFEL_vfns_fact_to_ren(self): ), }, ] - for fact_to_ren_lin, apfel_vals_dict in zip( - fact_to_ren_lin_list, apfel_vals_dict_list + for fact_to_ren_lin, apfel_vals_dict, nf_ref in zip( + fact_to_ren_lin_list, apfel_vals_dict_list, nf_refs ): # collect my values for order in [1, 2, 3]: From ba0a0cad21957dd112ee37c37adba2a50d59b486 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 1 Aug 2023 16:22:42 +0200 Subject: [PATCH 10/13] Update notebooks --- doc/source/overview/tutorials/alpha_s.ipynb | 12 +- doc/source/overview/tutorials/dglap.ipynb | 15 +-- doc/source/overview/tutorials/output.ipynb | 8 +- doc/source/overview/tutorials/pdf.ipynb | 130 ++++++++++---------- 4 files changed, 78 insertions(+), 87 deletions(-) diff --git a/doc/source/overview/tutorials/alpha_s.ipynb b/doc/source/overview/tutorials/alpha_s.ipynb index 327eb39c9..f9cc5a240 100644 --- a/doc/source/overview/tutorials/alpha_s.ipynb +++ b/doc/source/overview/tutorials/alpha_s.ipynb @@ -39,12 +39,10 @@ "from eko.quantities.heavy_quarks import MatchingScales, QuarkMassScheme\n", "\n", "# set the (alpha_s, alpha_em) reference values\n", - "couplings_ref = CouplingsInfo(\n", - " alphas=0.118, alphaem=0.007496252, ref=(91.0, 5)\n", - ")\n", + "couplings_ref = CouplingsInfo(alphas=0.118, alphaem=0.007496252, ref=(91.0, 5))\n", "\n", "# set heavy quark masses and their threshold ratios\n", - "heavy_quark_masses = np.power([1.51, 4.92, 172.0],2)\n", + "heavy_quark_masses = np.power([1.51, 4.92, 172.0], 2)\n", "thresholds_ratios = np.array([1.0, 1.0, 1.0])\n", "\n", "# set (QCD,QED) perturbative order\n", @@ -88,9 +86,9 @@ } ], "source": [ - "target_scale = 10.0 ** 2\n", + "target_scale = 10.0**2\n", "a_s = sc.a_s(target_scale)\n", - "print(\"The value of alpha_s at Q^2=100 GeV^2 is: \", 4. * np.pi * a_s)" + "print(\"The value of alpha_s at Q^2=100 GeV^2 is: \", 4.0 * np.pi * a_s)" ] }, { @@ -106,7 +104,7 @@ ], "metadata": { "kernelspec": { - "display_name": "eko-KkPVjVhh-py3.10", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, diff --git a/doc/source/overview/tutorials/dglap.ipynb b/doc/source/overview/tutorials/dglap.ipynb index 0e271a00c..98cb4a6e1 100644 --- a/doc/source/overview/tutorials/dglap.ipynb +++ b/doc/source/overview/tutorials/dglap.ipynb @@ -68,7 +68,7 @@ "th_card = example.theory()\n", "op_card = example.operator()\n", "# here we replace the grid with a very minimal one, to speed up the example\n", - "op_card.xgrid = [1e-3, 1e-2, 1e-1, 5e-1, 1.]" + "op_card.xgrid = [1e-3, 1e-2, 1e-1, 5e-1, 1.0]" ] }, { @@ -91,14 +91,9 @@ "{'order': [1, 0],\n", " 'couplings': {'alphas': 0.118,\n", " 'alphaem': 0.007496252,\n", - " 'scale': 91.2,\n", - " 'max_num_flavs': 6,\n", - " 'num_flavs_ref': 5,\n", + " 'ref': [91.2, 5],\n", " 'em_running': False},\n", - " 'heavy': {'num_flavs_init': 4,\n", - " 'num_flavs_max_pdf': 6,\n", - " 'intrinsic_flavors': [4],\n", - " 'masses': [[2.0, nan], [4.5, nan], [173.07, nan]],\n", + " 'heavy': {'masses': [[2.0, nan], [4.5, nan], [173.07, nan]],\n", " 'masses_scheme': 'pole',\n", " 'matching_ratios': [1.0, 1.0, 1.0]},\n", " 'xif': 1.0,\n", @@ -123,7 +118,7 @@ { "data": { "text/plain": [ - "{'mu0': 1.65,\n", + "{'init': [1.65, 4],\n", " 'mugrid': [[100.0, 5]],\n", " 'xgrid': [0.001, 0.01, 0.1, 0.5, 1.0],\n", " 'configs': {'evolution_method': 'iterate-exact',\n", @@ -174,7 +169,7 @@ "id": "880aadcf-8f87-4918-a0bc-09581d0d3579", "metadata": {}, "source": [ - "The actual result is a complicate EKO object, we will discuss it in a separate tutorial.\n", + "The actual result is a complicate EKO object, which we will discuss it in a separate tutorial.\n", "\n", "You have just run your first DGLAP calculation!" ] diff --git a/doc/source/overview/tutorials/output.ipynb b/doc/source/overview/tutorials/output.ipynb index 47345ac35..1cfe9f071 100644 --- a/doc/source/overview/tutorials/output.ipynb +++ b/doc/source/overview/tutorials/output.ipynb @@ -34,7 +34,7 @@ "id": "2f8f0666-c6ab-40f6-86f5-15773f205b51", "metadata": {}, "source": [ - "We can access the operator, by using the `open` method (similar to python's `open`):" + "We can access the operator, by using the `read` method:" ] }, { @@ -76,8 +76,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "TheoryCard(order=(1, 0), couplings=CouplingsInfo(alphas=0.118, alphaem=0.007496252, scale=91.2, max_num_flavs=6, num_flavs_ref=5, em_running=False), heavy=HeavyInfo(num_flavs_init=4, num_flavs_max_pdf=6, intrinsic_flavors=[4, 5, 6], masses=[[2.0, nan], [4.5, nan], [173.07, nan]], masses_scheme=, matching_ratios=[1.0, 1.0, 1.0]), xif=1.0, n3lo_ad_variation=(0, 0, 0, 0))\n", - "OperatorCard(mu0=1.65, mugrid=[(100.0, 5)], xgrid=, configs=Configs(evolution_method=, ev_op_max_order=(10, 0), ev_op_iterations=10, scvar_method=None, inversion_method=None, interpolation_polynomial_degree=4, interpolation_is_log=True, polarized=False, time_like=False, n_integration_cores=0), debug=Debug(skip_singlet=False, skip_non_singlet=False), eko_version='0.0.0')\n" + "TheoryCard(order=(1, 0), couplings=CouplingsInfo(alphas=0.118, alphaem=0.007496252, ref=(91.2, 5), em_running=False), heavy=HeavyInfo(masses=[[2.0, nan], [4.5, nan], [173.07, nan]], masses_scheme=, matching_ratios=[1.0, 1.0, 1.0]), xif=1.0, n3lo_ad_variation=(0, 0, 0, 0))\n", + "OperatorCard(init=(1.65, 4), mugrid=[(100.0, 5)], xgrid=, configs=Configs(evolution_method=, ev_op_max_order=(10, 0), ev_op_iterations=10, scvar_method=None, inversion_method=None, interpolation_polynomial_degree=4, interpolation_is_log=True, polarized=False, time_like=False, n_integration_cores=0), debug=Debug(skip_singlet=False, skip_non_singlet=False), eko_version='0.0.0')\n" ] } ], @@ -143,7 +143,7 @@ ], "source": [ "with eko.EKO.read(\"./myeko.tar\") as evolution_operator:\n", - " with evolution_operator.operator((10000.,5)) as op:\n", + " with evolution_operator.operator((10000.0, 5)) as op:\n", " print(f\"operator: {op.operator.shape}\")\n", " print(f\"error: {op.error.shape}\")" ] diff --git a/doc/source/overview/tutorials/pdf.ipynb b/doc/source/overview/tutorials/pdf.ipynb index 923520f59..bfc6ad3ad 100644 --- a/doc/source/overview/tutorials/pdf.ipynb +++ b/doc/source/overview/tutorials/pdf.ipynb @@ -13,9 +13,9 @@ "id": "3e0dcd0f", "metadata": {}, "source": [ - "## Method 1: Using apply_pdf\n", + "## Method 1: Using `apply_pdf`\n", "\n", - "In this first part we will compute the eko and subsequently apply the initial PDF \"manually\" calling a dedicated function. " + "In this first part, we compute the eko and subsequently apply the initial PDF \"manually\" calling a dedicated function. " ] }, { @@ -52,6 +52,7 @@ "import pathlib\n", "import eko\n", "from banana import toy\n", + "\n", "pdf = toy.mkPDF(\"\", 0)" ] }, @@ -71,6 +72,7 @@ "outputs": [], "source": [ "from ekobox.apply import apply_pdf\n", + "\n", "with eko.EKO.read(\"./myeko.tar\") as evolution_operator:\n", " evolved_pdfs = apply_pdf(evolution_operator, pdf)" ] @@ -92,7 +94,7 @@ { "data": { "text/plain": [ - "dict_keys([10000.0])" + "dict_keys([(10000.0, 5)])" ] }, "execution_count": 3, @@ -134,7 +136,7 @@ } ], "source": [ - "evolved_pdfs[10000.0][\"pdfs\"][21]" + "evolved_pdfs[(10000.0, 5)][\"pdfs\"][21]" ] }, { @@ -150,7 +152,7 @@ "id": "e925d2c9", "metadata": {}, "source": [ - "## Method 2: Using evolve_pdfs\n", + "## Method 2: Using `evolve_pdfs`\n", "\n", "In this second part we illustrate how to get (and install) directly a LHAPDF set evolved with eko. " ] @@ -171,6 +173,7 @@ "outputs": [], "source": [ "from banana import toy\n", + "\n", "pdf = toy.mkPDF(\"\", 0)" ] }, @@ -195,10 +198,10 @@ "th_card = example.theory()\n", "op_card = example.operator()\n", "# here we replace the grid with a very minimal one, to speed up the example\n", - "op_card.xgrid = eko.interpolation.XGrid([1e-3, 1e-2, 1e-1, 5e-1, 1.])\n", - "op_card.mugrid = [(10.,5), (100.,5)]\n", + "op_card.xgrid = eko.interpolation.XGrid([1e-3, 1e-2, 1e-1, 5e-1, 1.0])\n", + "op_card.mugrid = [(10.0, 5), (100.0, 5)]\n", "# set QCD LO evolution\n", - "th_card.orders = (1,0)" + "th_card.orders = (1, 0)" ] }, { @@ -236,16 +239,10 @@ ], "source": [ "from ekobox.evol_pdf import evolve_pdfs\n", + "\n", "path = pathlib.Path(\"./myeko2.tar\")\n", "path.unlink(missing_ok=True)\n", - "evolve_pdfs(\n", - " [pdf],\n", - " th_card,\n", - " op_card,\n", - " install=True,\n", - " name=\"Evolved_PDF\",\n", - " store_path=path\n", - ")" + "evolve_pdfs([pdf], th_card, op_card, install=True, name=\"Evolved_PDF\", store_path=path)" ] }, { @@ -273,6 +270,7 @@ ], "source": [ "import lhapdf\n", + "\n", "evolved_pdf = lhapdf.mkPDF(\"Evolved_PDF\", 0)" ] }, @@ -300,12 +298,12 @@ } ], "source": [ - "pid = 21 # gluon pid\n", - "Q2 = 89.10 # Q^2 in Gev^2\n", - "x = 0.01 # momentum fraction \n", + "pid = 21 # gluon pid\n", + "Q2 = 89.10 # Q^2 in Gev^2\n", + "x = 0.01 # momentum fraction\n", "\n", "# check that the particle is present\n", - "print(\"has gluon?\",evolved_pdf.hasFlavor(pid))\n", + "print(\"has gluon?\", evolved_pdf.hasFlavor(pid))\n", "# now do the lookup\n", "xg = evolved_pdf.xfxQ2(pid, x, Q2)\n", "print(f\"xg(x={x}, Q2={Q2}) = {xg}\")" @@ -352,24 +350,28 @@ "import lhapdf\n", "from ekobox.cards import example\n", "from eko.interpolation import make_grid\n", - "from eko.quantities.heavy_quarks import QuarkMassRef,HeavyQuarks\n", + "from eko.quantities.heavy_quarks import QuarkMassRef, HeavyQuarks\n", "\n", "# get the PDF object\n", "ct14llo = lhapdf.mkPDF(\"CT14llo\")\n", "\n", "# setup the operator card\n", "op_card = example.operator()\n", - "op_card.xgrid = eko.interpolation.XGrid(make_grid(30, 30)) # x grid\n", - "op_card.mugrid = [(float(q),5) for q in np.geomspace(5., 100, 5)] # Q2 grid\n", - "op_card.init = (1.295000,3) # starting point for the evolution \n", + "op_card.xgrid = eko.interpolation.XGrid(make_grid(30, 30)) # x grid\n", + "op_card.mugrid = [(float(q), 5) for q in np.geomspace(5.0, 100, 5)] # Q2 grid\n", + "op_card.init = (1.295000, 3) # starting point for the evolution\n", "\n", "# setup the theory card - this can be mostly inferred from the PDF's .info file\n", - "\n", "th_card = example.theory()\n", - "th_card.orders = (1,0) # QCD LO\n", - "th_card.heavy.masses = HeavyQuarks([QuarkMassRef([1.3,nan]), QuarkMassRef([4.75,nan]), QuarkMassRef([172.,nan])]) # quark mass\n", - "th_card.couplings.alphas = 0.130000 # reference value of alpha_s\n", - "th_card.couplings.ref = (91.1876,5) # the reference scale together with the number of flavors at which alpha_s is provided" + "th_card.orders = (1, 0) # QCD LO\n", + "th_card.heavy.masses = HeavyQuarks(\n", + " [QuarkMassRef([1.3, nan]), QuarkMassRef([4.75, nan]), QuarkMassRef([172.0, nan])]\n", + ") # quark mass\n", + "th_card.couplings.alphas = 0.130000 # reference value of alpha_s\n", + "th_card.couplings.ref = (\n", + " 91.1876,\n", + " 5,\n", + ") # the reference scale together with the number of flavors at which alpha_s is provided" ] }, { @@ -403,15 +405,11 @@ ], "source": [ "from ekobox.evol_pdf import evolve_pdfs\n", + "\n", "path = pathlib.Path(\"./myeko_ct14llo.tar\")\n", "path.unlink(missing_ok=True)\n", "evolve_pdfs(\n", - " [ct14llo],\n", - " th_card,\n", - " op_card,\n", - " install=True,\n", - " name=\"my_ct14llo\",\n", - " store_path=path\n", + " [ct14llo], th_card, op_card, install=True, name=\"my_ct14llo\", store_path=path\n", ")" ] }, @@ -435,33 +433,33 @@ "output_type": "stream", "text": [ "LHAPDF 6.4.0 loading /home/felix/local/share/LHAPDF/my_ct14llo/my_ct14llo_0000.dat\n", + "my_ct14llo PDF set, member #0, version 1\n", " x Q2 ct14llo my_ct14llo relative_diff\n", - "0 0.000010 25.000000 7.635785e+01 7.630461e+01 0.000697\n", - "1 0.000173 25.000000 3.194273e+01 3.192092e+01 0.000683\n", - "2 0.003000 25.000000 1.081843e+01 1.081086e+01 0.000701\n", - "3 0.051962 25.000000 1.958956e+00 1.958632e+00 0.000166\n", - "4 0.900000 25.000000 1.922415e-05 1.955026e-05 -0.016963\n", - "5 0.000010 111.803399 1.333957e+02 1.332985e+02 0.000729\n", - "6 0.000173 111.803399 4.777286e+01 4.773664e+01 0.000758\n", - "7 0.003000 111.803399 1.341028e+01 1.339967e+01 0.000791\n", - "8 0.051962 111.803399 1.978216e+00 1.978130e+00 0.000044\n", - "9 0.900000 111.803399 6.644805e-06 6.753652e-06 -0.016381\n", - "10 0.000010 500.000000 1.967032e+02 1.965456e+02 0.000801\n", - "11 0.000173 500.000000 6.291393e+01 6.286095e+01 0.000842\n", - "12 0.003000 500.000000 1.542347e+01 1.540996e+01 0.000876\n", - "13 0.051962 500.000000 1.947465e+00 1.947391e+00 0.000038\n", - "14 0.900000 500.000000 2.929060e-06 2.977306e-06 -0.016471\n", - "15 0.000010 2236.067977 2.633266e+02 2.631109e+02 0.000819\n", - "16 0.000173 2236.067977 7.708540e+01 7.701938e+01 0.000856\n", - "17 0.003000 2236.067977 1.700410e+01 1.698928e+01 0.000872\n", - "18 0.051962 2236.067977 1.893923e+00 1.893971e+00 -0.000025\n", - "19 0.900000 2236.067977 1.544450e-06 1.570997e-06 -0.017189\n", - "20 0.000010 10000.000000 3.314097e+02 3.311351e+02 0.000829\n", - "21 0.000173 10000.000000 9.023010e+01 9.015279e+01 0.000857\n", - "22 0.003000 10000.000000 1.825934e+01 1.824402e+01 0.000839\n", - "23 0.051962 10000.000000 1.830992e+00 1.831183e+00 -0.000104\n", - "24 0.900000 10000.000000 9.288458e-07 9.447927e-07 -0.017169\n", - "my_ct14llo PDF set, member #0, version 1\n" + "0 0.000010 25.000000 7.635785e+01 7.630719e+01 0.000663\n", + "1 0.000173 25.000000 3.194273e+01 3.192239e+01 0.000637\n", + "2 0.003000 25.000000 1.081843e+01 1.081160e+01 0.000632\n", + "3 0.051962 25.000000 1.958956e+00 1.958820e+00 0.000069\n", + "4 0.900000 25.000000 1.922415e-05 1.955440e-05 -0.017179\n", + "5 0.000010 111.803399 1.333957e+02 1.333028e+02 0.000697\n", + "6 0.000173 111.803399 4.777286e+01 4.773855e+01 0.000718\n", + "7 0.003000 111.803399 1.341028e+01 1.340044e+01 0.000734\n", + "8 0.051962 111.803399 1.978216e+00 1.978292e+00 -0.000038\n", + "9 0.900000 111.803399 6.644805e-06 6.756354e-06 -0.016787\n", + "10 0.000010 500.000000 1.967032e+02 1.965517e+02 0.000770\n", + "11 0.000173 500.000000 6.291393e+01 6.286327e+01 0.000805\n", + "12 0.003000 500.000000 1.542347e+01 1.541073e+01 0.000826\n", + "13 0.051962 500.000000 1.947465e+00 1.947532e+00 -0.000034\n", + "14 0.900000 500.000000 2.929060e-06 2.979511e-06 -0.017224\n", + "15 0.000010 2236.067977 2.633266e+02 2.631189e+02 0.000789\n", + "16 0.000173 2236.067977 7.708540e+01 7.702204e+01 0.000822\n", + "17 0.003000 2236.067977 1.700410e+01 1.699004e+01 0.000827\n", + "18 0.051962 2236.067977 1.893923e+00 1.894094e+00 -0.000090\n", + "19 0.900000 2236.067977 1.544450e-06 1.572860e-06 -0.018395\n", + "20 0.000010 10000.000000 3.314097e+02 3.311450e+02 0.000799\n", + "21 0.000173 10000.000000 9.023010e+01 9.015576e+01 0.000824\n", + "22 0.003000 10000.000000 1.825934e+01 1.824477e+01 0.000798\n", + "23 0.051962 10000.000000 1.830992e+00 1.831291e+00 -0.000163\n", + "24 0.900000 10000.000000 9.288458e-07 9.463689e-07 -0.018866\n" ] } ], @@ -471,15 +469,15 @@ "# load evolved pdf\n", "my_ct14llo = lhapdf.mkPDF(\"my_ct14llo\", 0)\n", "\n", - "pid = 21 # gluon pid\n", + "pid = 21 # gluon pid\n", "\n", "# collect data\n", - "log = {\"x\": [], \"Q2\" : [], \"ct14llo\": [], \"my_ct14llo\": [], \"relative_diff\": []} \n", - "for q in np.geomspace(5., 100, 5):\n", - " q2 = q**2.\n", + "log = {\"x\": [], \"Q2\": [], \"ct14llo\": [], \"my_ct14llo\": [], \"relative_diff\": []}\n", + "for q in np.geomspace(5.0, 100, 5):\n", + " q2 = q**2.0\n", " for x in np.geomspace(1e-5, 0.9, 5):\n", " value = ct14llo.xfxQ2(pid, x, q2)\n", - " my_value = my_ct14llo.xfxQ2(pid, x, q2)\n", + " my_value = my_ct14llo.xfxQ2(pid, x, q2)\n", " log[\"x\"].append(x)\n", " log[\"Q2\"].append(q2)\n", " log[\"ct14llo\"].append(value)\n", From 7cf9e1b30d2f17553094e23db15775aabc8e9545 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 1 Aug 2023 17:57:50 +0200 Subject: [PATCH 11/13] Use correct types in LHA --- benchmarks/lha_paper_bench.py | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/benchmarks/lha_paper_bench.py b/benchmarks/lha_paper_bench.py index 9254307dc..f92755af9 100644 --- a/benchmarks/lha_paper_bench.py +++ b/benchmarks/lha_paper_bench.py @@ -12,19 +12,15 @@ register(__file__) +_sqrt2 = float(np.sqrt(2.0)) + base_theory = { "ModEv": "EXA", - "Q0": np.sqrt( - 2.0 - ), # Eq. (30) :cite:`Giele:2002hx`, Eq. (4.53) :cite:`Dittmar:2005ed` - "mc": np.sqrt( - 2.0 - ), # Eq. (34) :cite:`Giele:2002hx`, Eq. (4.56) :cite:`Dittmar:2005ed` + "Q0": _sqrt2, # Eq. (30) :cite:`Giele:2002hx`, Eq. (4.53) :cite:`Dittmar:2005ed` + "mc": _sqrt2, # Eq. (34) :cite:`Giele:2002hx`, Eq. (4.56) :cite:`Dittmar:2005ed` "mb": 4.5, "mt": 175, - "Qref": np.sqrt( - 2.0 - ), # Eq. (32) :cite:`Giele:2002hx`,Eq. (4.53) :cite:`Dittmar:2005ed` + "Qref": _sqrt2, # Eq. (32) :cite:`Giele:2002hx`,Eq. (4.53) :cite:`Dittmar:2005ed` "alphas": 0.35, # Eq. (4.55) :cite:`Dittmar:2005ed` "alphaqed": 0.007496, "QED": 0, @@ -78,11 +74,11 @@ def sv_theories(self, pto): """ low = self.theory.copy() low["PTO"] = pto - low["XIF"] = np.sqrt(1.0 / 2.0) + low["XIF"] = 1.0 / _sqrt2 low["ModSV"] = "exponentiated" high = self.theory.copy() high["PTO"] = pto - high["XIF"] = np.sqrt(2.0) + high["XIF"] = _sqrt2 high["ModSV"] = "exponentiated" return [high, low] @@ -301,7 +297,7 @@ def benchmark_sv(self, pto): sv_theory["kcThr"] = 1.0 + 1e-15 sv_theory["nfref"] = 4 sv_theory["EScaleVar"] = 0 - low["XIR"] = np.sqrt(2.0) - high["XIR"] = np.sqrt(0.5) + low["XIR"] = _sqrt2 + high["XIR"] = 1.0 / _sqrt2 self.run_lha([low, high]) From 7bf06e2ffaa0cdb892004e29cb838d40b9786bbc Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 15 Jul 2024 16:08:38 +0300 Subject: [PATCH 12/13] Fix outdated sc tests --- tests/eko/scale_variations/test_diff.py | 4 +--- tests/eko/scale_variations/test_expanded.py | 4 +--- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/tests/eko/scale_variations/test_diff.py b/tests/eko/scale_variations/test_diff.py index c0dd67236..6c169e3bd 100644 --- a/tests/eko/scale_variations/test_diff.py +++ b/tests/eko/scale_variations/test_diff.py @@ -27,9 +27,7 @@ def compute_a_s(q2, order): couplings=CouplingsInfo( alphas=0.1181, alphaem=0.007496, - scale=91.00, - max_num_flavs=4, - num_flavs_ref=4, + ref=(91.00, 4), ), order=order, method=CouplingEvolutionMethod.EXPANDED, diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 6af22dd51..be14e45a6 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -18,9 +18,7 @@ def compute_a_s(q2, order): couplings=CouplingsInfo( alphas=0.1181, alphaem=0.007496, - scale=91.00, - max_num_flavs=4, - num_flavs_ref=4, + ref=(91.00, 4), ), order=order, method=CouplingEvolutionMethod.EXPANDED, From 618639b606b1a7aa93276dcc1fd9e87c51c51ba1 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 15 Jul 2024 16:12:20 +0300 Subject: [PATCH 13/13] Appease black --- src/ekomark/data/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ekomark/data/__init__.py b/src/ekomark/data/__init__.py index b6242845b..89fb6377b 100644 --- a/src/ekomark/data/__init__.py +++ b/src/ekomark/data/__init__.py @@ -1,4 +1,5 @@ """EKO database configuration.""" + from typing import Union from eko.io.runcards import Legacy, OperatorCard, TheoryCard