diff --git a/benchmarks/apfel_bench.py b/benchmarks/apfel_bench.py index 882b982ac..68f087233 100644 --- a/benchmarks/apfel_bench.py +++ b/benchmarks/apfel_bench.py @@ -38,14 +38,14 @@ def skip_pdfs(_theory): class BenchmarkVFNS(ApfelBenchmark): - """Benckmark VFNS""" + """Benchmark VFNS""" vfns_theory = { "FNS": "ZM-VFNS", "ModEv": [ "EXA", - # "EXP", - # "TRN", + "EXP", + "TRN", ], "kcThr": 1.0, "kbThr": 1.0, @@ -53,6 +53,9 @@ class BenchmarkVFNS(ApfelBenchmark): "Qref": np.sqrt(2.0), "alphas": 0.35, "Q0": np.sqrt(2.0), + "nfref": 3, + "nf0": 3, + "mc": 1.51, } vfns_theory = tolist(vfns_theory) @@ -65,11 +68,20 @@ def benchmark_plain(self, pto): cartesian_product(th), operators.build(operators.apfel_config), ["ToyLH"] ) - def benchmark_sv(self, pto): + def benchmark_sv(self, pto, svmode): """Scale Variation""" th = self.vfns_theory.copy() - th.update({"PTO": [pto], "XIR": [0.7071067811865475, 1.4142135623730951]}) + th.update( + { + "PTO": [pto], + "XIR": [1 / np.sqrt(2.0)], + "fact_to_ren_scale_ratio": [np.sqrt(2.0)], + "ModSV": [svmode], + "EScaleVar": [0], + "nfref": [4], + } + ) self.run( cartesian_product(th), operators.build(operators.apfel_config), ["ToyLH"] ) @@ -147,7 +159,7 @@ def benchmark_plain(self, pto): cartesian_product(th), operators.build(operators.apfel_config), ["ToyLH"] ) - def benchmark_sv(self, pto): + def benchmark_sv(self, pto, svmode): """Scale Variation""" ts = [] @@ -157,6 +169,8 @@ def benchmark_sv(self, pto): "PTO": [pto], "XIR": [np.sqrt(0.5)], "fact_to_ren_scale_ratio": [np.sqrt(2.0)], + "ModSV": [svmode], + "EScaleVar": [0], } ) ts.extend(cartesian_product(th)) @@ -166,6 +180,8 @@ def benchmark_sv(self, pto): "PTO": [pto], "XIR": [np.sqrt(2.0)], "fact_to_ren_scale_ratio": [np.sqrt(0.5)], + "ModSV": [svmode], + "EScaleVar": [0], } ) ts.extend(cartesian_product(th)) @@ -177,7 +193,7 @@ def benchmark_sv(self, pto): obj = BenchmarkVFNS() # obj = BenchmarkFFNS() - # obj.benchmark_plain(1) - # obj.benchmark_sv(1) + # obj.benchmark_plain(2) + obj.benchmark_sv(2, "exponentiated") # obj.benchmark_kthr(2) - obj.benchmark_msbar(2) + # obj.benchmark_msbar(2) diff --git a/benchmarks/banana.yaml b/benchmarks/banana.yaml index 98c209e7f..7c28b7cb9 100644 --- a/benchmarks/banana.yaml +++ b/benchmarks/banana.yaml @@ -1,5 +1,7 @@ -database_path: data/benchmark.db -input_tables: - - theories - - operators - - cache +paths: + database: data/benchmark.db +input: + tables: + - theories + - operators + - cache diff --git a/benchmarks/eko/benchmark_ad.py b/benchmarks/eko/benchmark_ad.py index 9b0ac352e..b220cae5a 100644 --- a/benchmarks/eko/benchmark_ad.py +++ b/benchmarks/eko/benchmark_ad.py @@ -14,7 +14,6 @@ def benchmark_melling_g3_pegasus(): check_melling_g3_pegasus(N) -@pytest.mark.isolated def check_melling_g3_pegasus(N): S1 = h.harmonic_S1(N) N1 = N + 1.0 @@ -40,25 +39,23 @@ def check_melling_g3_pegasus(N): - 0.3174 * (zeta2 - S15 / N5) / N5 + 0.0699 * (zeta2 - S16 / N6) / N6 ) - np.testing.assert_allclose(h.mellin_g3(N), SPMOM) @pytest.mark.isolated def benchmark_gamma_ns_1_pegasus(): - for N in range(10): + # remember that singlet has pole at N=1 + for N in [2, 3, 4, 1 + 1j, 1 - 1j, 2 + 1j, 3 + 1j]: for NF in [3, 4, 5]: - check_gamma_ns_1_pegasus(N, NF) + check_gamma_1_pegasus(N, NF) -@pytest.mark.isolated -def check_gamma_ns_1_pegasus(N, NF): - # pylint: disable=line-too-long,too-many-locals,too-many-statements +def check_gamma_1_pegasus(N, NF): # Test against pegasus implementation ZETA2 = h.zeta2 ZETA3 = h.zeta3 - N = np.random.rand(1) + np.random.rand(1) * 1j + # N = np.random.rand(1) + np.random.rand(1) * 1j S1 = h.harmonic_S1(N) S2 = h.harmonic_S2(N) diff --git a/benchmarks/eko/benchmark_evol_to_unity.py b/benchmarks/eko/benchmark_evol_to_unity.py index bd427a018..df5540b34 100644 --- a/benchmarks/eko/benchmark_evol_to_unity.py +++ b/benchmarks/eko/benchmark_evol_to_unity.py @@ -34,6 +34,7 @@ class BenchmarkBackwardForward: "MaxNfPdf": 6, "MaxNfAs": 6, "HQ": "POLE", + "ModSV": None, } operators_card = { "Q2grid": [10], diff --git a/benchmarks/eko/benchmark_msbar_evolution.py b/benchmarks/eko/benchmark_msbar_evolution.py index 73e94c325..da8db85e4 100644 --- a/benchmarks/eko/benchmark_msbar_evolution.py +++ b/benchmarks/eko/benchmark_msbar_evolution.py @@ -1,9 +1,12 @@ # -*- coding: utf-8 -*- """This module benchmarks MSbar mass evolution against APFEL.""" +import copy + import numpy as np import pytest from eko import msbar_masses +from eko.basis_rotation import quark_names from eko.strong_coupling import StrongCoupling # try to load APFEL - if not available, we'll use the cached values @@ -40,12 +43,18 @@ @pytest.mark.isolated class BenchmarkMSbar: def benchmark_APFEL_msbar_evolution(self): - Q2s = np.power([1, 96, 150], 2) - alphas_ref = 0.118 - scale_ref = 91.0**2 - thresholds_ratios = np.power((1.0, 1.0, 1.0), 2) - Q2m = np.power([2.0, 4.5, 175], 2) - m2 = np.power((1.4, 4.5, 175), 2) + theory = copy.deepcopy(theory_dict) + bench_values = dict(zip(np.power([1, 96, 150], 2), [3, 5, 5])) + theory.update( + { + "mc": 1.4, + "mb": 4.5, + "mt": 175.0, + "Qmc": 2.0, + "Qmb": 4.5, + "Qmt": 175.0, + } + ) apfel_vals_dict = { "exact": { 0: np.array( @@ -96,27 +105,32 @@ def benchmark_APFEL_msbar_evolution(self): } # collect my values for method in ["exact", "expanded"]: + theory["ModEv"] = "EXP" if method == "expanded" else "EXA" for order in [0, 1, 2]: + theory["PTO"] = order as_VFNS = StrongCoupling( - alphas_ref, - scale_ref, - m2, - thresholds_ratios, + theory["alphas"], + theory["Qref"] ** 2, + msbar_masses.compute(theory), + np.power([theory["kcThr"], theory["kbThr"], theory["ktThr"]], 2), order=order, method=method, - hqm_scheme="MSBAR", + nf_ref=theory["nfref"], + hqm_scheme=theory["HQ"], ) my_vals = [] - for Q2 in Q2s: + for Q2, nf_to in bench_values.items(): my_masses = [] for n in [3, 4, 5]: my_masses.append( msbar_masses.evolve( - m2[n - 3], - Q2m[n - 3], + theory[f"m{quark_names[n]}"] ** 2, + theory[f"Qm{quark_names[n]}"] ** 2, strong_coupling=as_VFNS, fact_to_ren=1.0, q2_to=Q2, + nf_ref=n, + nf_to=nf_to, ) ) my_vals.append(my_masses) @@ -128,15 +142,17 @@ def benchmark_APFEL_msbar_evolution(self): apfel.SetTheory("QCD") apfel.SetPerturbativeOrder(order) apfel.SetAlphaEvolution(method) - apfel.SetAlphaQCDRef(alphas_ref, np.sqrt(scale_ref)) + apfel.SetAlphaQCDRef(theory["alphas"], theory["Qref"]) apfel.SetVFNS() - apfel.SetMSbarMasses(*np.sqrt(m2)) - apfel.SetMassScaleReference(*np.sqrt(Q2m)) + apfel.SetMSbarMasses(theory["mc"], theory["mb"], theory["mt"]) + apfel.SetMassScaleReference( + theory["Qmc"], theory["Qmb"], theory["Qmt"] + ) apfel.SetRenFacRatio(1.0) apfel.InitializeAPFEL() # collect apfel masses apfel_vals_cur = [] - for Q2 in Q2s: + for Q2 in bench_values: masses = [] for n in [4, 5, 6]: masses.append(apfel.HeavyQuarkMass(n, np.sqrt(Q2))) @@ -147,7 +163,7 @@ def benchmark_APFEL_msbar_evolution(self): ) # check myself to APFEL np.testing.assert_allclose( - apfel_vals, np.sqrt(np.array(my_vals)), rtol=2e-3 + apfel_vals, np.sqrt(np.array(my_vals)), rtol=2.3e-3 ) def benchmark_APFEL_msbar_solution(self): @@ -184,10 +200,9 @@ def benchmark_APFEL_msbar_solution(self): ) apfel.SetRenFacRatio(theory_dict["fact_to_ren_scale_ratio"]) apfel.InitializeAPFEL() - apfel.EnableWelcomeMessage(1) # check myself to APFEL np.testing.assert_allclose( - apfel_vals, np.sqrt(np.array(my_masses)), rtol=5e-5 + apfel_vals, np.sqrt(np.array(my_masses)), rtol=4e-4 ) def benchmark_msbar_solution_kthr(self): diff --git a/benchmarks/eko/benchmark_strong_coupling.py b/benchmarks/eko/benchmark_strong_coupling.py index af57d3106..39bbbd517 100644 --- a/benchmarks/eko/benchmark_strong_coupling.py +++ b/benchmarks/eko/benchmark_strong_coupling.py @@ -787,3 +787,81 @@ def benchmark_lhapdf_zmvfns_lo(self): # Max absolute difference: 2.58611282e-06 # Max relative difference: 0.00013379 np.testing.assert_allclose(lhapdf_vals, np.array(my_vals), rtol=1.5e-4) + + def benchmark_APFEL_fact_to_ren_lha_settings(self): + + theory_dict = { + "alphas": 0.35, + "Qref": np.sqrt(2.0), + "nfref": 4, + "nf0": 3, + "MaxNfPdf": 6, + "MaxNfAs": 6, + "Q0": np.sqrt(2.0), + "fact_to_ren_scale_ratio": np.sqrt(2.0), + "mc": np.sqrt(2.0), + "mb": 4.5, + "mt": 175.0, + "kcThr": 1.0, + "kbThr": 1.0, + "ktThr": 1.0, + "HQ": "POLE", + "Qmc": np.sqrt(2.0), + "Qmb": 4.5, + "Qmt": 175.0, + "PTO": 2, + "ModEv": "EXA", + } + Q2s = [2.0] + sc = StrongCoupling.from_dict(theory_dict) + fact_to_ren = theory_dict["fact_to_ren_scale_ratio"] ** 2 + for Q2 in Q2s: + + my_val = sc.a_s(Q2 / fact_to_ren, Q2) + path = sc.thresholds.path(Q2 / fact_to_ren) + my_val_4 = sc.a_s(Q2 / fact_to_ren, Q2, nf_to=4) + path_4 = sc.thresholds.path(Q2 / fact_to_ren, 4) + my_val_3 = sc.a_s(Q2 / fact_to_ren, Q2, nf_to=3) + path_3 = sc.thresholds.path(Q2 / fact_to_ren, 3) + + # path_4 it's not matched + assert len(path_4) == 1 + + # path_3 is the same as path backward in nf and in q2. + assert len(path_3) == 2 + assert len(path) == 2 + assert path_3[1].nf < path_3[0].nf + assert path_3[1].q2_from < path_3[0].q2_from + + apfel_val_ref = 0.03478112968976964 + if use_APFEL: + # run apfel + apfel.CleanUp() + apfel.SetTheory("QCD") + apfel.SetPerturbativeOrder(theory_dict["PTO"]) + apfel.SetAlphaEvolution("exact") + apfel.SetAlphaQCDRef(theory_dict["alphas"], theory_dict["Qref"]) + apfel.SetVFNS() + apfel.SetPoleMasses( + theory_dict["mc"], + theory_dict["mb"], + theory_dict["mt"], + ) + apfel.SetMassMatchingScales( + theory_dict["kcThr"], + theory_dict["kbThr"], + theory_dict["ktThr"], + ) + apfel.SetRenFacRatio(1.0 / theory_dict["fact_to_ren_scale_ratio"]) + apfel.InitializeAPFEL() + # collect a_s + apfel_val = apfel.AlphaQCD( + np.sqrt(Q2) / theory_dict["fact_to_ren_scale_ratio"] + ) / (4.0 * np.pi) + # check APFEL cached value + np.testing.assert_allclose(apfel_val_ref, apfel_val) + + # check myself to APFEL + np.testing.assert_allclose(apfel_val_ref, my_val, rtol=0.03) + np.testing.assert_allclose(apfel_val_ref, my_val_4) + np.testing.assert_allclose(my_val, my_val_3) diff --git a/benchmarks/lha_paper_bench.py b/benchmarks/lha_paper_bench.py index 9ef572229..415ed665d 100644 --- a/benchmarks/lha_paper_bench.py +++ b/benchmarks/lha_paper_bench.py @@ -3,15 +3,18 @@ Benchmark to :cite:`Giele:2002hx` (LO + NLO) and :cite:`Dittmar:2005ed` (NNLO) """ import numpy as np +from banana import register +from eko.interpolation import make_lambert_grid from ekomark.benchmark.runner import Runner +register(__file__) + base_theory = { "ModEv": "EXA", "Q0": np.sqrt( 2.0 ), # Eq. (30) :cite:`Giele:2002hx`, Eq. (4.53) :cite:`Dittmar:2005ed` - "nfref": 3, "mc": np.sqrt( 2.0 ), # Eq. (34) :cite:`Giele:2002hx`, Eq. (4.56) :cite:`Dittmar:2005ed` @@ -34,11 +37,11 @@ class LHABenchmark(Runner): Globally set the external program to LHA """ - external = "LHA" - - theory = {} - - rotate_to_evolution_basis = True + def __init__(self): + super().__init__() + self.external = "LHA" + self.theory = {} + self.rotate_to_evolution_basis = True def plain_theory(self, pto): """ @@ -75,9 +78,11 @@ def sv_theories(self, pto): low = self.theory.copy() low["PTO"] = pto low["fact_to_ren_scale_ratio"] = np.sqrt(1.0 / 2.0) + low["ModSV"] = "exponentiated" high = self.theory.copy() high["PTO"] = pto high["fact_to_ren_scale_ratio"] = np.sqrt(2.0) + high["ModSV"] = "exponentiated" return [high, low] @staticmethod @@ -112,7 +117,7 @@ def run_lha(self, theory_updates): { "Q2grid": [1e4], "ev_op_iterations": 10, - # "debug_skip_singlet": True + "interpolation_xgrid": make_lambert_grid(60).tolist(), } ], ["ToyLH"], @@ -130,32 +135,37 @@ def benchmark_sv(self, pto): class BenchmarkVFNS(LHABenchmark): """Variable Flavor Number Scheme""" - theory = base_theory.copy() - theory.update( - { - "FNS": "ZM-VFNS", # ignored by eko, but needed by LHA_utils - "kcThr": 1.0, - "kbThr": 1.0, - "ktThr": 1.0, - "nf0": 3, - } - ) + def __init__(self): + super().__init__() + self.theory = base_theory.copy() + self.theory.update( + { + "FNS": "ZM-VFNS", # ignored by eko, but needed by LHA_utils + "kcThr": 1.0, + "kbThr": 1.0, + "ktThr": 1.0, + "nf0": 3, + "nfref": 3, + } + ) class BenchmarkFFNS(LHABenchmark): """Fixed Flavor Number Scheme""" - theory = base_theory.copy() - theory.update( - { - "FNS": "FFNS", # ignored by eko, but needed by LHA_utils - "NfFF": 4, - "nfref": 4, - "kcThr": 0.0, - "kbThr": np.inf, - "ktThr": np.inf, - } - ) + def __init__(self): + super().__init__() + self.theory = base_theory.copy() + self.theory.update( + { + "FNS": "FFNS", # ignored by eko, but needed by LHA_utils + "NfFF": 4, + "nfref": 4, + "kcThr": 0.0, + "kbThr": np.inf, + "ktThr": np.inf, + } + ) @staticmethod def skip_pdfs(theory): @@ -168,10 +178,58 @@ def skip_pdfs(theory): return ffns_skip_pdfs +class BenchmarkRunner(BenchmarkVFNS): + """ + Generic benchmark runner using the LHA VFNS settings + """ + + def __init__(self, external): + super().__init__() + self.external = external + self.sandbox = True + + def benchmark_sv(self, pto): + """ + Scale variations + + Parameters + ---------- + pto : int + perturbation order + + Returns + ------- + list(dict) + theory updates + """ + high, low = self.sv_theories(pto) + + # here we need to adjust also the + # apfel initial nf, which can't + # be accessed in other ways + if self.external == "apfel": + for sv_theory in [low, high]: + 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) + + self.run_lha([low, high]) + + if __name__ == "__main__": - obj = BenchmarkVFNS() + # Benchmark to LHA + # obj = BenchmarkVFNS() # obj = BenchmarkFFNS() - obj.benchmark_plain(2) - # obj.benchmark_sv(1) + # obj.benchmark_plain(2) + # obj.benchmark_sv(2) + + # # VFNS benchmarks with LHA settings + programs = ["LHA", "pegasus", "apfel"] + for p in programs: + obj = BenchmarkRunner(p) + # obj.benchmark_plain(2) + obj.benchmark_sv(2) diff --git a/benchmarks/pegasus_bench.py b/benchmarks/pegasus_bench.py index 36fcc3a9f..498ba322f 100644 --- a/benchmarks/pegasus_bench.py +++ b/benchmarks/pegasus_bench.py @@ -3,11 +3,14 @@ Benchmark to Pegasus :cite:`Vogt:2004ns` """ import numpy as np +from banana import register from banana.data import cartesian_product from ekomark.benchmark.runner import Runner from ekomark.data import operators +register(__file__) + def tolist(input_dict): output_dict = input_dict.copy() @@ -54,6 +57,8 @@ class BenchmarkVFNS(PegasusBenchmark): "Qref": np.sqrt(2.0), "alphas": 0.35, "Q0": np.sqrt(2.0), + "nfref": 3, + "nf0": 3, } zm_theory = tolist(zm_theory) @@ -70,7 +75,7 @@ def benchmark_plain(self, pto): cartesian_product(th), operators.build(operators.pegasus_config), ["ToyLH"] ) - def benchmark_sv(self, pto): + def benchmark_sv(self, pto, svmode): """Scale Variation""" th = self.zm_theory.copy() @@ -78,6 +83,7 @@ def benchmark_sv(self, pto): { "PTO": [pto], "fact_to_ren_scale_ratio": [np.sqrt(0.5), np.sqrt(2.0)], + "ModSV": [svmode], } ) self.run( @@ -118,7 +124,7 @@ def benchmark_plain(self, pto): cartesian_product(th), operators.build(operators.pegasus_config), ["ToyLH"] ) - def benchmark_sv(self, pto): + def benchmark_sv(self, pto, svmode): """Scale Variation""" th = self.ffns_theory.copy() @@ -126,6 +132,7 @@ def benchmark_sv(self, pto): { "PTO": [pto], "fact_to_ren_scale_ratio": [np.sqrt(0.5), np.sqrt(2.0)], + "ModSV": [svmode], } ) self.run( @@ -137,7 +144,7 @@ def benchmark_sv(self, pto): obj = BenchmarkVFNS() # obj = BenchmarkFFNS() - obj.benchmark_plain(1) + # obj.benchmark_plain(1) - # obj.benchmark_sv(1) + obj.benchmark_sv(2, "exponentiated") # vfns.benchmark_sv() diff --git a/benchmarks/sandbox.py b/benchmarks/sandbox.py index 9ca7d6f6a..4cb511cbd 100644 --- a/benchmarks/sandbox.py +++ b/benchmarks/sandbox.py @@ -1,8 +1,8 @@ # -*- coding: utf-8 -*- # pylint: skip-file import numpy as np +from banana import register -from ekomark import register from ekomark.benchmark.runner import Runner from ekomark.data import operators diff --git a/doc/.gitignore b/doc/.gitignore index 3025cff24..d60bb3d69 100644 --- a/doc/.gitignore +++ b/doc/.gitignore @@ -3,6 +3,7 @@ source/development/code_todos.rst # Ignore auto generated module references source/modules source/development/ekomark +source/code/ekobox # ignore temporary build files _build/ # Ignore generated sphinx-bibtex file diff --git a/doc/Makefile b/doc/Makefile index 57aadc8d8..3954fb760 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -10,7 +10,9 @@ BUILDDIR = build EKODIR = ../src/eko EKOOUT = $(SOURCEDIR)/modules/eko EKOMARKDIR = ../src/ekomark -EKOMARKOUT = ./source/development/ekomark +EKOMARKOUT = $(SOURCEDIR)/development/ekomark +EKOBOXKDIR = ../src/ekobox +EKOBOXKOUT = $(SOURCEDIR)/code/ekobox TODOOUTFILE = ./source/development/code_todos.rst @@ -35,7 +37,9 @@ clean: rm -rf _build cleanall: clean clean-todos - rm -rf $(SOURCEDIR)/modules + rm -rf $(EKOOUT) + rm -rf $(EKOMARKOUT) + rm -rf $(EKOBOXKOUT) # TODOs todos: diff --git a/doc/source/code/IO-tabs/InterpolatorDispatcher.rst b/doc/source/code/IO-tabs/InterpolatorDispatcher.rst deleted file mode 100644 index f8c128cc7..000000000 --- a/doc/source/code/IO-tabs/InterpolatorDispatcher.rst +++ /dev/null @@ -1,5 +0,0 @@ -InterpolatorDispatcher parameter: - -- :py:obj:`list(float)` ``interpolation_xgrid`` the interpolation grid; **Required** -- :py:obj:`int` ``interpolation_polynomial_degree`` polynomial degree of the interpolating function; **Required** -- :py:obj:`bool` ``interpolation_is_log`` use logarithmic interpolation? **Required** diff --git a/doc/source/code/IO-tabs/ThresholdConfig.rst b/doc/source/code/IO-tabs/ThresholdConfig.rst deleted file mode 100644 index cc1f4393c..000000000 --- a/doc/source/code/IO-tabs/ThresholdConfig.rst +++ /dev/null @@ -1,6 +0,0 @@ -ThresholdConfig parameters: - -- :py:obj:`float` ``Q0`` reference scale from which to start; **Required** -- :py:obj:`float` ``mc`` charm mass in GeV; **Required** -- :py:obj:`float` ``mb`` bottom mass in GeV; **Required** -- :py:obj:`float` ``mt`` top mass in GeV; **Required** diff --git a/doc/source/code/Operators.rst b/doc/source/code/Operators.rst index 80f48003f..c58a64b37 100644 --- a/doc/source/code/Operators.rst +++ b/doc/source/code/Operators.rst @@ -55,7 +55,7 @@ The classes are nested as follows: * is the connection of the :class:`~eko.evolution_operator.Operator` between the different flavor bases * is initialized with the 3-dimensional :ref:`theory/FlavorSpace:Operator Anomalous Dimension Basis` - * does recombine the operator in the :ref:`theory/FlavorSpace:Operator Intrinsic Evolution Basis` + * does recombine the operator in the :ref:`theory/FlavorSpace:Operator Intrinsic QCD Evolution Basis` (see :doc:`Matching Conditions `) * exports the operators to :ref:`theory/FlavorSpace:Operator Flavor Basis` in a :class:`~numpy.ndarray` @@ -63,4 +63,4 @@ The classes are nested as follows: * represents a single operator in Mellin space for a given element of the :ref:`theory/FlavorSpace:Operator Bases` * inside :class:`~eko.evolution_operator.Operator` they are in :ref:`theory/FlavorSpace:Operator Anomalous Dimension Basis` - * inside :class:`~eko.evolution_operator.physical.PhysicalOperator` they are in :ref:`theory/FlavorSpace:Operator Intrinsic Evolution Basis` + * inside :class:`~eko.evolution_operator.physical.PhysicalOperator` they are in :ref:`theory/FlavorSpace:Operator Intrinsic QCD Evolution Basis` diff --git a/doc/source/code/Utilities.rst b/doc/source/code/Utilities.rst index d19a9262f..2c7a6b83e 100644 --- a/doc/source/code/Utilities.rst +++ b/doc/source/code/Utilities.rst @@ -8,10 +8,69 @@ Apart from the :doc:`operator ` classes, `eko` also provides some uti - Implementation of the flavor number scheme and the quark thresholds both for the :class:`eko.strong_coupling.StrongCoupling` and the :doc:`operators <../theory/Matching>` + When running in |VFNS| it is important to specify the number of flavors active at each given scale, since the evolution path + can be different depending of the chosen setting. This path is determined by :meth:`eko.thresholds.ThresholdsAtlas.path`. -.. include:: IO-tabs/ThresholdConfig.rst + Let us consider two examples to better illustrate how it works. + Imagine to have a boundary condition ``q2_ref=1``, ``nf_ref=3`` with heavy quarks mass thresholds + at: ``mc=2``, ``mb=3``, ``mt=4`` and would like to evolve your object (|PDF| or :math:`a_s`) to an higher + scale (say ``q2_to=49``). The corresponding ``ThresholdsAtlas`` will look like: + .. code-block:: -- :class:`eko.strong_coupling.StrongCoupling` Implementation of the :ref:`theory/pQCD:strong coupling` + ThresholdsAtlas [0.00e+00 - 4.00e+00 - 9.00e+00 - 16.00e+00 - inf], ref=1.000000000000001 @ 3 -- :class:`eko.interpolation.InterpolatorDispatcher` Implementation of the :doc:`../theory/Interpolation` + and the evolution path will be given by: + + .. code-block:: + + [PathSegment(1 -> 4, nf=3), PathSegment(4 -> 9, nf=4), PathSegment(9 -> 16, nf=5), PathSegment(16 -> 49, nf=6)] + + where automatically the number of flavors active in each :class:`eko.thresholds.PathSegment` is determined according with the + specified thresholds scales. + + However some more complicated situations can arise when the boundary conditions are not given with the same prescription + of the :class:`eko.thresholds.ThresholdAtlas`. Let's now consider as boundary condition ``q2_ref=64`` with ``nf_ref=3``: + + .. code-block:: + + ThresholdsAtlas [0.00e+00 - 4.00e+00 - 9.00e+00 - 16.00e+00 - inf], ref=64.000000000000001 @ 3 + + Again we would like to evolve to ``q2_to=49`` but now giving a different number of active flavors to the final scale. + Some possible paths are: + + - a path to ``q2_to=49`` with ``nf=6`` which is given by the following list of :class:`eko.thresholds.PathSegment` : + + .. code-block:: + + [PathSegment(64 -> 4, nf=3), PathSegment(4 -> 9, nf=4), PathSegment(9 -> 16, nf=5), PathSegment(16 -> 49, nf=6)] + + and it's determined according to the list of heavy quark thresholds. This path is default option when the argument + ``nf_to`` is not given. + + - A path to ``q2_to=49`` enforcing ``nf=4``: + + .. code-block:: + + [PathSegment(64 -> 4, nf=3), PathSegment(4 -> 49, nf=4)] + + - A path to ``q2_to=49`` enforcing ``nf=3`` which will simply contain a single :class:`eko.thresholds.PathSegment`: + + .. code-block:: + + [PathSegment(64 -> 49, nf=3)] + + + In the first two cases as first step, you go back to the closest matching scale (closest in `nf`), + running in ``nf=3``, since that is what has been specified in the boundary conditions. + Then you are back in the flavor flow, and you cross all the other thresholds going according to the prescription given + by the :class:`eko.thresholds.ThresholdAtlas`. + While in the third example you go directly to the final scale, since there is no matching scale in the middle. + You can notice that :meth:`eko.thresholds.ThresholdsAtlas.path` are always ordered according to `nf` and not `q2` scales. + This property is used in |VFNS| to determine if the matchings are done for an increasing or decreasing number of + light flavors. + + +- :class:`eko.strong_coupling.StrongCoupling`: Implementation of the :ref:`theory/pQCD:strong coupling` + +- :class:`eko.interpolation.InterpolatorDispatcher`: Implementation of the :doc:`../theory/Interpolation` diff --git a/doc/source/code/genpdf.rst b/doc/source/code/genpdf.rst new file mode 100644 index 000000000..551dee766 --- /dev/null +++ b/doc/source/code/genpdf.rst @@ -0,0 +1,75 @@ +genpdf +====== + +We provide also a console script called ``genpdf`` that is able +to generate and install a custom |PDF| set in the `lhapdf` format. +In particular, the command ``genpdf install [NAME]`` simply +install the |PDF| called ``[NAME]`` in the lhapdf folder and +the ``genpdf generate [NAME]`` command generates the custom |PDF| +and saves it as ``[NAME]`` in the current directory. +Notice that the argument ``[NAME]`` is the only mandatory one. + +The custom |PDF| can be generated in three different ways which +are accessible trough the option ``-p [PARENT]`` (the complete spelling +is ``genpdf generate [NAME] -p [PARENT]``): + + 1. If ``[PARENT]`` is the name of an available |PDF|, it is used as parent + |PDF| and thus copied to generate the new custom PDF. + 2. If ``[PARENT]`` is "toylh" or "toy", the **toy** |PDF| is used as parent. + 3. If the option ``-p [PARENT]`` is not used, the |PDF| is + generated using **x(1-x)** for all the flavors. + +Trough the use of the argument +``[LABEL]`` (``genpdf generate [NAME] [LABEL] [OPTIONS]``) it is also possible +to specify a set of flavors (expressed in |pid| basis) or a set of +**evolution basis components** on which filtering the custom |PDF|. +In this way the specified set is kept in the final |PDF|, while the rest +is discarded. + +In the case of custom |PDF| generated starting from a parent |PDF|, +it is possible to generate all the members trough the flag ``-m``. If this +flag is not used, only the *zero* member is generated (together with the *info* +file of course). Using the flag ``-m`` when the custom |PDF| is generated +either using the **toy** |PDF| or the **x(1-x)** function, has no effects. + +In order to automatically install the custom |PDF| in the lhapdf folder +at generation time (so without using ``genpdf install [NAME]`` after the +generation), it is possible to use the ``-i`` flag. + + +We also provide an API with some additional features and possibilities +such as generating a |PDF| with a custom function for every |pid| +(through a ``dict`` structure) and filtering custom combination of +flavors - see :mod:`here ` for details. + +Examples +-------- + +.. code-block:: bash + + $ genpdf generate gonly 21 + +This will generate the custom PDF using the debug x(1-x) PDF as parent +and then it will keep only the gluon. + +.. code-block:: bash + + $ genpdf install gonly + +This will install the previous PDF in the lhapdf folder. + +.. code-block:: bash + + $ genpdf generate Sonly S -p toy -i + +This will generate the custom PDF using the toy PDF as parent and then +it will keep only the singlet combination. The generated PDF is also +automatically installed in the lhapdf folder. + +.. code-block:: bash + + $ genpdf generate Vonly V -p CT10 -m + +This will generate the custom PDF using the CT10 PDF set as parent +(if available) and it will keep only the valence combination. Moreover +it will generate all the members of the parent PDF. diff --git a/doc/source/code/interface.rst b/doc/source/code/interface.rst new file mode 100644 index 000000000..18f9c42b6 --- /dev/null +++ b/doc/source/code/interface.rst @@ -0,0 +1,9 @@ +Interface +========= + + +.. toctree:: + :maxdepth: 1 + + ekobox/ekobox.rst + genpdf.rst diff --git a/doc/source/conf.py b/doc/source/conf.py index efc04dbf3..706b6937b 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -305,10 +305,9 @@ def process_numba_docstring( """Recover the docstring under numba, as the numba.njit decorator doesn't repeat the __doc__""" if not isinstance(obj, nb.core.registry.CPUDispatcher): return - else: - original = obj.py_func - orig_sig = inspect.signature(original) - lines = orig_sig.__doc__ + original = obj.py_func + orig_sig = inspect.signature(original) + lines = orig_sig.__doc__ # https://github.com/readthedocs/readthedocs.org/issues/1139#issuecomment-312626491 @@ -319,15 +318,14 @@ def run_apidoc(_): sys.path.append(str(here.parent)) # 'eko' - docs_dest = here / "modules" / "eko" - package = here.parents[1] / "src" / "eko" - main(["--module-first", "-o", str(docs_dest), str(package)]) - (docs_dest / "modules.rst").unlink() - # 'ekomark' - docs_dest = here / "development" / "ekomark" - package = here.parents[1] / "src" / "ekomark" - main(["--module-first", "-o", str(docs_dest), str(package)]) - (docs_dest / "modules.rst").unlink() + for pkg, docs_dest in dict( + eko=here / "modules" / "eko", + ekomark=here / "development" / "ekomark", + ekobox=here / "code" / "ekobox", + ).items(): + package = here.parents[1] / "src" / pkg + main(["--module-first", "-o", str(docs_dest), str(package)]) + (docs_dest / "modules.rst").unlink() def setup(app): diff --git a/doc/source/development/ekomark.rst b/doc/source/development/ekomark.rst index 9a9d99b4d..5a8630645 100644 --- a/doc/source/development/ekomark.rst +++ b/doc/source/development/ekomark.rst @@ -30,7 +30,7 @@ Ekomark is composed by four subpackages: * ``plot`` containing all the scripts to produce the output plots. -The banana configuration is loaded from ``banana_cfg.py`` file. +The banana configuration is loaded from ``banana/cfg.py`` file. To run Ekomark see the section of the available :doc:`runners`. Furthermore Ekomark provides also a python interpreter called `ekonavigator` to inspect the cached benchmark results. diff --git a/doc/source/index.rst b/doc/source/index.rst index f0c78d36f..e724727d5 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -46,6 +46,7 @@ EKO is ... :hidden: code/IO + code/interface API code/Operators code/Utilities diff --git a/doc/source/shared/abbreviations.rst b/doc/source/shared/abbreviations.rst index e7742731e..017ce001f 100644 --- a/doc/source/shared/abbreviations.rst +++ b/doc/source/shared/abbreviations.rst @@ -40,14 +40,19 @@ .. |pid| replace:: :abbr:`PID ((Monte Carlo) parton identifier)` +.. QCD + .. |OME| replace:: - :abbr:`OME (operator matrix element)` + :abbr:`OME (Operator Matrix Element)` .. |MSbar| replace:: :math:`\overline{MS}` .. |RGE| replace:: - :abbr:`RGE (renormalization group equation)` + :abbr:`RGE (Renormalization Group Equation)` + +.. |MHOU| replace:: + :abbr:`MHOU (Missing Higher Order Uncertainties)` .. |QCD| replace:: :abbr:`QCD (Quantum Chromodynamics)` @@ -89,3 +94,7 @@ .. |T| raw:: html ✓ + + +.. |API| replace:: + :abbr:`API (Application Program Interface)` diff --git a/doc/source/theory/Matching.rst b/doc/source/theory/Matching.rst index b2e2046a1..cc2902022 100644 --- a/doc/source/theory/Matching.rst +++ b/doc/source/theory/Matching.rst @@ -66,8 +66,7 @@ Furthermore in the right side basis :math:`\tilde h^{+}, \tilde h^{-}` are intri The :math:`\mathbf{A}^{(n_f)}(\mu_{h+1}^2)` can be computed order by order in :math:`a_s`: .. math :: - \mathbf{A}^{(n_f)}(\mu_{h}^2) = \mathbf{I} + a_s^{(n_f)}(\mu_{h}^2) \mathbf{A}^{(n_f),(1)} + \left(a_s^{(n_f)}(\mu_{h}^2)\right)^2 \mathbf{A}^{(n_f),(2)} - + \mathbf{A}^{(n_f)}(\mu_{h}^2) = \mathbf{I} + a_s^{(n_f+1)}(\mu_{h}^2) \mathbf{A}^{(n_f),(1)} + \left(a_s^{(n_f+1)}(\mu_{h}^2)\right)^2 \mathbf{A}^{(n_f),(2)} where the :math:`\mathbf{A}^{(n_f),(k)}` are given up to |NNLO| by the following expressions: @@ -83,6 +82,8 @@ The coefficients :math:`A^{(n_f),(k)}_{ij}(z,\mu_{h}^2)` have been firstly compu been :doc:`Mellin transformed ` to be used inside EKO. They depend on the scale :math:`\mu_{h}^2` only through the logarithm :math:`\ln(\mu_{h}^2/m_{h}^2)`, in particular the coefficient :math:`A_{gg,H}^{S,(1)}` is fully proportional to :math:`\ln(\mu_{h}^2/m_{h}^2)`. +During the matching we use :math:`a_s^{(n_f+1)}`: in fact the :math:`a_s` decoupling gives raise to some additional logarithms +:math:`\ln(\mu_{h}^2/m_{h}^2)`, which are cancelled by the OME's :math:`A_{kl,H}`. We remark that contributions of the higher quark at |NNLO| have not been computed yet, thus the elements :math:`A_{qH}^{(2)},A_{gH}^{(2)}A_{HH}^{(2)}` are not encoded in EKO despite of being present. On the other hand the elements :math:`A_{qq}^{ps},A_{qg}` are known to start at |N3LO|. diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst index b101c6280..2611c9876 100644 --- a/doc/source/theory/pQCD.rst +++ b/doc/source/theory/pQCD.rst @@ -46,6 +46,11 @@ have to be applied :cite:`Schroder:2005hy,Chetyrkin:2005ia`. In particular, the matching involved in the change from :math:`n_f` to :math:`n_f-1` schemes is presented in equation 3.1 of :cite:`Schroder:2005hy` for |MSbar| masses, while the same expression for POLE masses is reported in Appendix A. +For this reason the boundary conditions of :class:`eko.strong_coupling.StrongCoupling` +can be specified at ``scale_ref`` along with ``nf_ref`` and, the computed result can +depend on the number of flavors at the target scale, see :meth:`eko.strong_coupling.StrongCoupling.a_s` +An example how the evolution path is determined is given :doc:`here`. + QCD Splitting Functions ----------------------- @@ -80,23 +85,85 @@ The expression of the pure |QED| and of the mixed |QED| :math:`\otimes` |QCD| sp Scale Variations ---------------- -The usual procedure in solving |DGLAP| that is also applied :doc:`here +The usual procedure in solving |DGLAP| applied :doc:`here ` is to rewrite the equations in term of the running coupling :math:`a_s` assuming the factorization scale :math:`\mu_F^2` (the inherit scale of the |PDF|) and the renormalization scale :math:`\mu_R^2` (the inherit scale -for the strong coupling) to be equal. This constraint can however be lifted by a -suitable redefinition of the splitting kernels :cite:`Vogt:2004ns`: +for the strong coupling) to be equal. +This constraint, however, can be lifted in order to provide an estimation of the +missing higher order uncertainties (|MHOU|) coming from |DGLAP| evolution :cite:`AbdulKhalek:2019ihb`. +Since scale-dependent contributions to a perturbative prediction are fixed by |RGE| invariance, +the scale variation can be used to generate higher order contributions, +which are then taken as a proxy for the whole missing higher orders. +This method provides many advantages: -.. math :: - \gamma^{(1)}(N) &\to \gamma^{(1)}(N) - \beta_0 \ln(\mu_F^2/\mu_R^2) \gamma^{(0)} \\ - \gamma^{(2)}(N) &\to \gamma^{(2)}(N) - 2 \beta_0 \ln(\mu_F^2/\mu_R^2) \gamma^{(1)} - ( \beta_1 \ln(\mu_F^2/\mu_R^2) - \beta_0^2 \ln^2(\mu_F^2/\mu_R^2) ) \gamma^{(0)} + * it naturally incorporates renormalization group invariance, + as the perturbative order increases, estimates of |MHOU| decrease; + * the same procedure can be used for any perturbative process, + since the scale dependence of the strong coupling :math:`a_s(\mu^2)` and of PDFs is universal; + +However, there is no unique prescription to determine the specific range of the scale variation, +the most used prescription specify to vary the factor :math:`\mu_F/\mu_R` in the range: +:math:`1/2 \le \mu_F/\mu_R \le 2`. + +This variation can be performed at least at two different levels during the |PDF| +evolution, always evaluating the strong coupling at :math:`\mu_R^2`. + + * For ``ModSV='exponentiated'`` the variation is applied directly to the splitting functions + and the anomalous dimension are then modified using :cite:`Vogt:2004ns`: + + .. math :: + & \gamma^{(1)}(N) \to \gamma^{(1)}(N) - \beta_0 k \gamma^{(0)} \\ + & \gamma^{(2)}(N) \to \gamma^{(2)}(N) - 2 \beta_0 k \gamma^{(1)} - ( \beta_1 k - \beta_0^2 k^2) \gamma^{(0)} \\ + & \gamma^{(3)}(N) \to \gamma^{(3)}(N) - 3 \beta_0 k \gamma^{(2)} - ( 2 \beta_1 k - 3 \beta_0^2 k^2) \gamma^{(1)} - (\beta_2 k - \frac{5}{2} \beta_1 \beta_0 k^2 + \beta_0^3 k^3) \gamma^{(0)} \\ + & k = \ln(\mu_F^2/\mu_R^2) + + This procedure corresponds to Eq. (3.32) of :cite:`AbdulKhalek:2019ihb`, and we recommend to use it along with + ``ModEv='iterate-exact'`` in order to be in agreement with the treatment of the evolution integral expansion. + + + * In ``ModSV='expanded'`` the |EKO| is multiplied by an additional kernel, such that + the scale variation is applied to the whole |PDF| set: + + .. math :: + \tilde{\mathbf{E}}(a_s \leftarrow a_s^0) & = \tilde{\mathbf{K}}(a_s) \tilde{\mathbf{E}}(a_s \leftarrow a_s^0) \\ + \tilde{\mathbf{K}}(a_s) & = 1 - k \gamma + \frac{1}{2} k^2 \left ( \gamma^{2} - \beta \frac{\partial \gamma}{\partial a_s} \right ) \\ + & + \frac{1}{6} k^3 \left [ - \beta \frac{\partial}{\partial a_s} \left( \beta \frac{\partial \gamma}{\partial a_s} \right) + 3 \beta \frac{\partial \gamma}{\partial a_s} \gamma - \gamma^3 \right ] + \mathcal{O}(k^4) + + where scale variation kernel is expanded consistently order by order in :math:`a_s`, + leading to: + + .. math :: + \tilde{\mathbf{K}}(a_s) \approx & 1 - a_s k \gamma^{(0)} + a_s^2 \left [ - k \gamma^{(1)} + \frac{1}{2} k^2 \gamma^{(0)} (\beta_0 + \gamma^{(0)}) \right ] \\ + & + a_s^3 \left [ -k \gamma^{(2)} + \frac{1}{2} k^2 \left(\beta_1 \gamma^{(0)} + 2 \gamma^{(1)} (\beta_0 + \gamma^{(0)} ) \right) \right. \\ + & \left. - \frac{1}{6} k^3 \gamma^{(0)} \left(2 \beta_0^2 + 3 \beta_0 \gamma^{(0)}+\gamma^{(0),2} \right) \right] + \mathcal{O}(a^4) + + + In this way the dependence of the |EKO| on :math:`k` is factorized outside the unvaried evolution kernel. + This procedure is repeated for each different flavor patch present in the evolution path. + It corresponds to Eq. (3.35) of :cite:`AbdulKhalek:2019ihb`, and we recommend to use it along with + ``ModEv='truncated'`` in order to keep consistency with the evolution integral expansion. + + +By construction, the corrections of the order :math:`\mathcal{O}(k^n)` will appear +at the order :math:`n` in the expansion :math:`a_s`. +This happens because :math:`\beta \approx \mathcal{O}(a_s^2)`, :math:`\gamma \approx \mathcal{O}(a_s)` +and the contribution proportional to :math:`\mathcal{O}(k^n)` is originated +by the `n-th` derivative in :math:`\gamma` :cite:`AbdulKhalek:2019ihb`. +Furthermore the distance between the varied |EKO| and the unvaried one will decrease while +keeping higher order terms in :math:`a_s` -while keeping the evaluation of the strong coupling always at :math:`\mu_R^2`. -Estimating the theoretical uncertainties imposed on |PDF| determination due to -missing higher order corrections using scale variation in the evolution -corresponds to schemes A and B in :cite:`AbdulKhalek:2019ihb`. +Notice that in principle the two methods should be equivalent, especially for fully +linearized solutions (``ModEv=truncated``, ``ev_op_iterations=1``), +where the difference depends only on the perturbative expansion in :math:`a_s`. +However, in our implementation this is not exactly true; +since the integral of :math:`-\frac{\gamma(a_s)}{\beta(a_s)}` is evaluated before +the scale variation procedure is applied, the difference between the two schemes +depends also on the actual evolution distance and on the ratio :math:`k`. +When using the scale variations, boundary conditions for the strong coupling :math:`a_s` +(``Qref`` and ``nfref``) have to be given according to renormalization scales. Heavy Quark Masses ------------------ diff --git a/poetry.lock b/poetry.lock index 708b1e231..4594bf08a 100644 --- a/poetry.lock +++ b/poetry.lock @@ -122,28 +122,6 @@ SQLAlchemy = ">=1.4.29,<2.0.0" [package.extras] docs = ["Sphinx (>=4.3.2,<5.0.0)", "sphinx-rtd-theme (>=1.0.0,<2.0.0)", "sphinxcontrib-bibtex (>=2.4.1,<3.0.0)"] -[[package]] -name = "black" -version = "22.1.0" -description = "The uncompromising code formatter." -category = "dev" -optional = false -python-versions = ">=3.6.2" - -[package.dependencies] -click = ">=8.0.0" -mypy-extensions = ">=0.4.3" -pathspec = ">=0.9.0" -platformdirs = ">=2" -tomli = ">=1.1.0" -typing-extensions = {version = ">=3.10.0.0", markers = "python_version < \"3.10\""} - -[package.extras] -colorama = ["colorama (>=0.4.3)"] -d = ["aiohttp (>=3.7.4)"] -jupyter = ["ipython (>=7.8.0)", "tokenize-rt (>=3.2.0)"] -uvloop = ["uvloop (>=0.15.2)"] - [[package]] name = "certifi" version = "2021.10.8" @@ -176,7 +154,7 @@ name = "click" version = "8.0.4" description = "Composable command line interface toolkit" category = "main" -optional = false +optional = true python-versions = ">=3.6" [package.dependencies] @@ -534,14 +512,6 @@ category = "dev" optional = false python-versions = "*" -[[package]] -name = "mypy-extensions" -version = "0.4.3" -description = "Experimental type system extensions for programs checked with the mypy typechecker." -category = "dev" -optional = false -python-versions = "*" - [[package]] name = "nodeenv" version = "1.6.0" @@ -614,14 +584,6 @@ python-versions = ">=3.6" qa = ["flake8 (==3.8.3)", "mypy (==0.782)"] testing = ["docopt", "pytest (<6.0.0)"] -[[package]] -name = "pathspec" -version = "0.9.0" -description = "Utility library for gitignore style pattern matching of file paths." -category = "dev" -optional = false -python-versions = "!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*,!=3.4.*,>=2.7" - [[package]] name = "pdbpp" version = "0.10.3" @@ -1129,7 +1091,7 @@ test = ["pytest"] [[package]] name = "sqlalchemy" -version = "1.4.31" +version = "1.4.32" description = "Database Abstraction Library" category = "main" optional = true @@ -1225,7 +1187,7 @@ socks = ["PySocks (>=1.5.6,!=1.5.7,<2.0)"] [[package]] name = "virtualenv" -version = "20.13.2" +version = "20.13.3" description = "Virtual Python Environment builder" category = "dev" optional = false @@ -1285,7 +1247,7 @@ mark = ["banana-hep", "sqlalchemy", "pandas", "matplotlib"] [metadata] lock-version = "1.1" python-versions = "^3.8,<3.11" -content-hash = "5ba606f5ae173bbdf7ac42303b6cbcf36262d220aa8c3041ba6064bc3df26c0a" +content-hash = "1b6fb7851a767927d7a6b3aceafd6fcc9a7ee27752e5c402629b173fac56fea9" [metadata.files] a3b2bbc3ced97675ac3a71df45f55ba = [ @@ -1337,31 +1299,6 @@ banana-hep = [ {file = "banana-hep-0.6.3.tar.gz", hash = "sha256:296e5e7fa79c2271348258b6f2165cc34067b2a9c60e2c97aec035eb645fe014"}, {file = "banana_hep-0.6.3-py3-none-any.whl", hash = "sha256:aa52344c4fac9e9e5475dbe175da44fd52ab37bb25b68c0a71eff6f704cd5b82"}, ] -black = [ - {file = "black-22.1.0-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:1297c63b9e1b96a3d0da2d85d11cd9bf8664251fd69ddac068b98dc4f34f73b6"}, - {file = "black-22.1.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:2ff96450d3ad9ea499fc4c60e425a1439c2120cbbc1ab959ff20f7c76ec7e866"}, - {file = "black-22.1.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:0e21e1f1efa65a50e3960edd068b6ae6d64ad6235bd8bfea116a03b21836af71"}, - {file = "black-22.1.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:e2f69158a7d120fd641d1fa9a921d898e20d52e44a74a6fbbcc570a62a6bc8ab"}, - {file = "black-22.1.0-cp310-cp310-win_amd64.whl", hash = "sha256:228b5ae2c8e3d6227e4bde5920d2fc66cc3400fde7bcc74f480cb07ef0b570d5"}, - {file = "black-22.1.0-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:b1a5ed73ab4c482208d20434f700d514f66ffe2840f63a6252ecc43a9bc77e8a"}, - {file = "black-22.1.0-cp36-cp36m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:35944b7100af4a985abfcaa860b06af15590deb1f392f06c8683b4381e8eeaf0"}, - {file = "black-22.1.0-cp36-cp36m-win_amd64.whl", hash = "sha256:7835fee5238fc0a0baf6c9268fb816b5f5cd9b8793423a75e8cd663c48d073ba"}, - {file = "black-22.1.0-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:dae63f2dbf82882fa3b2a3c49c32bffe144970a573cd68d247af6560fc493ae1"}, - {file = "black-22.1.0-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:5fa1db02410b1924b6749c245ab38d30621564e658297484952f3d8a39fce7e8"}, - {file = "black-22.1.0-cp37-cp37m-win_amd64.whl", hash = "sha256:c8226f50b8c34a14608b848dc23a46e5d08397d009446353dad45e04af0c8e28"}, - {file = "black-22.1.0-cp38-cp38-macosx_10_9_universal2.whl", hash = "sha256:2d6f331c02f0f40aa51a22e479c8209d37fcd520c77721c034517d44eecf5912"}, - {file = "black-22.1.0-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:742ce9af3086e5bd07e58c8feb09dbb2b047b7f566eb5f5bc63fd455814979f3"}, - {file = "black-22.1.0-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:fdb8754b453fb15fad3f72cd9cad3e16776f0964d67cf30ebcbf10327a3777a3"}, - {file = "black-22.1.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:f5660feab44c2e3cb24b2419b998846cbb01c23c7fe645fee45087efa3da2d61"}, - {file = "black-22.1.0-cp38-cp38-win_amd64.whl", hash = "sha256:6f2f01381f91c1efb1451998bd65a129b3ed6f64f79663a55fe0e9b74a5f81fd"}, - {file = "black-22.1.0-cp39-cp39-macosx_10_9_universal2.whl", hash = "sha256:efbadd9b52c060a8fc3b9658744091cb33c31f830b3f074422ed27bad2b18e8f"}, - {file = "black-22.1.0-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:8871fcb4b447206904932b54b567923e5be802b9b19b744fdff092bd2f3118d0"}, - {file = "black-22.1.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:ccad888050f5393f0d6029deea2a33e5ae371fd182a697313bdbd835d3edaf9c"}, - {file = "black-22.1.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:07e5c049442d7ca1a2fc273c79d1aecbbf1bc858f62e8184abe1ad175c4f7cc2"}, - {file = "black-22.1.0-cp39-cp39-win_amd64.whl", hash = "sha256:373922fc66676133ddc3e754e4509196a8c392fec3f5ca4486673e685a421321"}, - {file = "black-22.1.0-py3-none-any.whl", hash = "sha256:3524739d76b6b3ed1132422bf9d82123cd1705086723bc3e235ca39fd21c667d"}, - {file = "black-22.1.0.tar.gz", hash = "sha256:a7c0192d35635f6fc1174be575cb7915e92e5dd629ee79fdaf0dcfa41a80afb5"}, -] certifi = [ {file = "certifi-2021.10.8-py2.py3-none-any.whl", hash = "sha256:d62a0163eb4c2344ac042ab2bdf75399a71a2d8c7d47eac2e2ee91b9d6339569"}, {file = "certifi-2021.10.8.tar.gz", hash = "sha256:78884e7c1d4b00ce3cea67b44566851c4343c120abd683433ce934a68ea58872"}, @@ -1775,10 +1712,6 @@ mccabe = [ {file = "mccabe-0.6.1-py2.py3-none-any.whl", hash = "sha256:ab8a6258860da4b6677da4bd2fe5dc2c659cff31b3ee4f7f5d64e79735b80d42"}, {file = "mccabe-0.6.1.tar.gz", hash = "sha256:dd8d182285a0fe56bace7f45b5e7d1a6ebcbf524e8f3bd87eb0f125271b8831f"}, ] -mypy-extensions = [ - {file = "mypy_extensions-0.4.3-py2.py3-none-any.whl", hash = "sha256:090fedd75945a69ae91ce1303b5824f428daf5a028d2f6ab8a299250a846f15d"}, - {file = "mypy_extensions-0.4.3.tar.gz", hash = "sha256:2d82818f5bb3e369420cb3c4060a7970edba416647068eb4c5343488a6c604a8"}, -] nodeenv = [ {file = "nodeenv-1.6.0-py2.py3-none-any.whl", hash = "sha256:621e6b7076565ddcacd2db0294c0381e01fd28945ab36bcf00f41c5daf63bef7"}, {file = "nodeenv-1.6.0.tar.gz", hash = "sha256:3ef13ff90291ba2a4a7a4ff9a979b63ffdd00a464dbe04acf0ea6471517a4c2b"}, @@ -1873,10 +1806,6 @@ parso = [ {file = "parso-0.8.3-py2.py3-none-any.whl", hash = "sha256:c001d4636cd3aecdaf33cbb40aebb59b094be2a74c556778ef5576c175e19e75"}, {file = "parso-0.8.3.tar.gz", hash = "sha256:8c07be290bb59f03588915921e29e8a50002acaf2cdc5fa0e0114f91709fafa0"}, ] -pathspec = [ - {file = "pathspec-0.9.0-py2.py3-none-any.whl", hash = "sha256:7d15c4ddb0b5c802d161efc417ec1a2558ea2653c2e8ad9c19098201dc1c993a"}, - {file = "pathspec-0.9.0.tar.gz", hash = "sha256:e564499435a2673d586f6b2130bb5b95f04a3ba06f81b8f895b651a3c76aabb1"}, -] pdbpp = [ {file = "pdbpp-0.10.3-py2.py3-none-any.whl", hash = "sha256:79580568e33eb3d6f6b462b1187f53e10cd8e4538f7d31495c9181e2cf9665d1"}, {file = "pdbpp-0.10.3.tar.gz", hash = "sha256:d9e43f4fda388eeb365f2887f4e7b66ac09dce9b6236b76f63616530e2f669f5"}, @@ -2145,42 +2074,41 @@ sphinxcontrib-serializinghtml = [ {file = "sphinxcontrib_serializinghtml-1.1.5-py2.py3-none-any.whl", hash = "sha256:352a9a00ae864471d3a7ead8d7d79f5fc0b57e8b3f95e9867eb9eb28999b92fd"}, ] sqlalchemy = [ - {file = "SQLAlchemy-1.4.31-cp27-cp27m-macosx_10_14_x86_64.whl", hash = "sha256:c3abc34fed19fdeaead0ced8cf56dd121f08198008c033596aa6aae7cc58f59f"}, - {file = "SQLAlchemy-1.4.31-cp27-cp27m-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:8d0949b11681380b4a50ac3cd075e4816afe9fa4a8c8ae006c1ca26f0fa40ad8"}, - {file = "SQLAlchemy-1.4.31-cp27-cp27m-win32.whl", hash = "sha256:f3b7ec97e68b68cb1f9ddb82eda17b418f19a034fa8380a0ac04e8fe01532875"}, - {file = "SQLAlchemy-1.4.31-cp27-cp27m-win_amd64.whl", hash = "sha256:81f2dd355b57770fdf292b54f3e0a9823ec27a543f947fa2eb4ec0df44f35f0d"}, - {file = "SQLAlchemy-1.4.31-cp27-cp27mu-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:4ad31cec8b49fd718470328ad9711f4dc703507d434fd45461096da0a7135ee0"}, - {file = "SQLAlchemy-1.4.31-cp310-cp310-macosx_10_15_x86_64.whl", hash = "sha256:05fa14f279d43df68964ad066f653193187909950aa0163320b728edfc400167"}, - {file = "SQLAlchemy-1.4.31-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:dccff41478050e823271642837b904d5f9bda3f5cf7d371ce163f00a694118d6"}, - {file = "SQLAlchemy-1.4.31-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:57205844f246bab9b666a32f59b046add8995c665d9ecb2b7b837b087df90639"}, - {file = "SQLAlchemy-1.4.31-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ea8210090a816d48a4291a47462bac750e3bc5c2442e6d64f7b8137a7c3f9ac5"}, - {file = "SQLAlchemy-1.4.31-cp310-cp310-win32.whl", hash = "sha256:2e216c13ecc7fcdcbb86bb3225425b3ed338e43a8810c7089ddb472676124b9b"}, - {file = "SQLAlchemy-1.4.31-cp310-cp310-win_amd64.whl", hash = "sha256:e3a86b59b6227ef72ffc10d4b23f0fe994bef64d4667eab4fb8cd43de4223bec"}, - {file = "SQLAlchemy-1.4.31-cp36-cp36m-macosx_10_14_x86_64.whl", hash = "sha256:2fd4d3ca64c41dae31228b80556ab55b6489275fb204827f6560b65f95692cf3"}, - {file = "SQLAlchemy-1.4.31-cp36-cp36m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:6f22c040d196f841168b1456e77c30a18a3dc16b336ddbc5a24ce01ab4e95ae0"}, - {file = "SQLAlchemy-1.4.31-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:c0c7171aa5a57e522a04a31b84798b6c926234cb559c0939840c3235cf068813"}, - {file = "SQLAlchemy-1.4.31-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d046a9aeba9bc53e88a41e58beb72b6205abb9a20f6c136161adf9128e589db5"}, - {file = "SQLAlchemy-1.4.31-cp36-cp36m-win32.whl", hash = "sha256:d86132922531f0dc5a4f424c7580a472a924dd737602638e704841c9cb24aea2"}, - {file = "SQLAlchemy-1.4.31-cp36-cp36m-win_amd64.whl", hash = "sha256:ca68c52e3cae491ace2bf39b35fef4ce26c192fd70b4cd90f040d419f70893b5"}, - {file = "SQLAlchemy-1.4.31-cp37-cp37m-macosx_10_14_x86_64.whl", hash = "sha256:cf2cd387409b12d0a8b801610d6336ee7d24043b6dd965950eaec09b73e7262f"}, - {file = "SQLAlchemy-1.4.31-cp37-cp37m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:bb4b15fb1f0aafa65cbdc62d3c2078bea1ceecbfccc9a1f23a2113c9ac1191fa"}, - {file = "SQLAlchemy-1.4.31-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:c317ddd7c586af350a6aef22b891e84b16bff1a27886ed5b30f15c1ed59caeaa"}, - {file = "SQLAlchemy-1.4.31-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:3c7ed6c69debaf6198fadb1c16ae1253a29a7670bbf0646f92582eb465a0b999"}, - {file = "SQLAlchemy-1.4.31-cp37-cp37m-win32.whl", hash = "sha256:6a01ec49ca54ce03bc14e10de55dfc64187a2194b3b0e5ac0fdbe9b24767e79e"}, - {file = "SQLAlchemy-1.4.31-cp37-cp37m-win_amd64.whl", hash = "sha256:330eb45395874cc7787214fdd4489e2afb931bc49e0a7a8f9cd56d6e9c5b1639"}, - {file = "SQLAlchemy-1.4.31-cp38-cp38-macosx_10_14_x86_64.whl", hash = "sha256:5e9c7b3567edbc2183607f7d9f3e7e89355b8f8984eec4d2cd1e1513c8f7b43f"}, - {file = "SQLAlchemy-1.4.31-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:de85c26a5a1c72e695ab0454e92f60213b4459b8d7c502e0be7a6369690eeb1a"}, - {file = "SQLAlchemy-1.4.31-cp38-cp38-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:975f5c0793892c634c4920057da0de3a48bbbbd0a5c86f5fcf2f2fedf41b76da"}, - {file = "SQLAlchemy-1.4.31-cp38-cp38-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d5c20c8415173b119762b6110af64448adccd4d11f273fb9f718a9865b88a99c"}, - {file = "SQLAlchemy-1.4.31-cp38-cp38-win32.whl", hash = "sha256:b35dca159c1c9fa8a5f9005e42133eed82705bf8e243da371a5e5826440e65ca"}, - {file = "SQLAlchemy-1.4.31-cp38-cp38-win_amd64.whl", hash = "sha256:b7b20c88873675903d6438d8b33fba027997193e274b9367421e610d9da76c08"}, - {file = "SQLAlchemy-1.4.31-cp39-cp39-macosx_10_15_x86_64.whl", hash = "sha256:85e4c244e1de056d48dae466e9baf9437980c19fcde493e0db1a0a986e6d75b4"}, - {file = "SQLAlchemy-1.4.31-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:e79e73d5ee24196d3057340e356e6254af4d10e1fc22d3207ea8342fc5ffb977"}, - {file = "SQLAlchemy-1.4.31-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:15a03261aa1e68f208e71ae3cd845b00063d242cbf8c87348a0c2c0fc6e1f2ac"}, - {file = "SQLAlchemy-1.4.31-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:0ddc5e5ccc0160e7ad190e5c61eb57560f38559e22586955f205e537cda26034"}, - {file = "SQLAlchemy-1.4.31-cp39-cp39-win32.whl", hash = "sha256:289465162b1fa1e7a982f8abe59d26a8331211cad4942e8031d2b7db1f75e649"}, - {file = "SQLAlchemy-1.4.31-cp39-cp39-win_amd64.whl", hash = "sha256:9e4fb2895b83993831ba2401b6404de953fdbfa9d7d4fa6a4756294a83bbc94f"}, - {file = "SQLAlchemy-1.4.31.tar.gz", hash = "sha256:582b59d1e5780a447aada22b461e50b404a9dc05768da1d87368ad8190468418"}, + {file = "SQLAlchemy-1.4.32-cp27-cp27m-macosx_10_14_x86_64.whl", hash = "sha256:4b2bcab3a914715d332ca783e9bda13bc570d8b9ef087563210ba63082c18c16"}, + {file = "SQLAlchemy-1.4.32-cp27-cp27m-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:159c2f69dd6efd28e894f261ffca1100690f28210f34cfcd70b895e0ea7a64f3"}, + {file = "SQLAlchemy-1.4.32-cp27-cp27m-win_amd64.whl", hash = "sha256:d7e483f4791fbda60e23926b098702340504f7684ce7e1fd2c1bf02029288423"}, + {file = "SQLAlchemy-1.4.32-cp27-cp27mu-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:4aa96e957141006181ca58e792e900ee511085b8dae06c2d08c00f108280fb8a"}, + {file = "SQLAlchemy-1.4.32-cp310-cp310-macosx_10_15_x86_64.whl", hash = "sha256:576684771456d02e24078047c2567025f2011977aa342063468577d94e194b00"}, + {file = "SQLAlchemy-1.4.32-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:fff677fa4522dafb5a5e2c0cf909790d5d367326321aeabc0dffc9047cb235bd"}, + {file = "SQLAlchemy-1.4.32-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:8679f9aba5ac22e7bce54ccd8a77641d3aea3e2d96e73e4356c887ebf8ff1082"}, + {file = "SQLAlchemy-1.4.32-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:c7046f7aa2db445daccc8424f50b47a66c4039c9f058246b43796aa818f8b751"}, + {file = "SQLAlchemy-1.4.32-cp310-cp310-win32.whl", hash = "sha256:bedd89c34ab62565d44745212814e4b57ef1c24ad4af9b29c504ce40f0dc6558"}, + {file = "SQLAlchemy-1.4.32-cp310-cp310-win_amd64.whl", hash = "sha256:199dc6d0068753b6a8c0bd3aceb86a3e782df118260ebc1fa981ea31ee054674"}, + {file = "SQLAlchemy-1.4.32-cp36-cp36m-macosx_10_14_x86_64.whl", hash = "sha256:8e1e5d96b744a4f91163290b01045430f3f32579e46d87282449e5b14d27d4ac"}, + {file = "SQLAlchemy-1.4.32-cp36-cp36m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:edfcf93fd92e2f9eef640b3a7a40db20fe3c1d7c2c74faa41424c63dead61b76"}, + {file = "SQLAlchemy-1.4.32-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:04164e0063feb7aedd9d073db0fd496edb244be40d46ea1f0d8990815e4b8c34"}, + {file = "SQLAlchemy-1.4.32-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:5ba59761c19b800bc2e1c9324da04d35ef51e4ee9621ff37534bc2290d258f71"}, + {file = "SQLAlchemy-1.4.32-cp36-cp36m-win32.whl", hash = "sha256:708973b5d9e1e441188124aaf13c121e5b03b6054c2df59b32219175a25aa13e"}, + {file = "SQLAlchemy-1.4.32-cp36-cp36m-win_amd64.whl", hash = "sha256:316270e5867566376e69a0ac738b863d41396e2b63274616817e1d34156dff0e"}, + {file = "SQLAlchemy-1.4.32-cp37-cp37m-macosx_10_14_x86_64.whl", hash = "sha256:9a0195af6b9050c9322a97cf07514f66fe511968e623ca87b2df5e3cf6349615"}, + {file = "SQLAlchemy-1.4.32-cp37-cp37m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:f7e4a3c0c3c596296b37f8427c467c8e4336dc8d50f8ed38042e8ba79507b2c9"}, + {file = "SQLAlchemy-1.4.32-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:bca714d831e5b8860c3ab134c93aec63d1a4f493bed20084f54e3ce9f0a3bf99"}, + {file = "SQLAlchemy-1.4.32-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:e9a680d9665f88346ed339888781f5236347933906c5a56348abb8261282ec48"}, + {file = "SQLAlchemy-1.4.32-cp37-cp37m-win32.whl", hash = "sha256:9cb5698c896fa72f88e7ef04ef62572faf56809093180771d9be8d9f2e264a13"}, + {file = "SQLAlchemy-1.4.32-cp37-cp37m-win_amd64.whl", hash = "sha256:8b9a395122770a6f08ebfd0321546d7379f43505882c7419d7886856a07caa13"}, + {file = "SQLAlchemy-1.4.32-cp38-cp38-macosx_10_14_x86_64.whl", hash = "sha256:3f88a4ee192142eeed3fe173f673ea6ab1f5a863810a9d85dbf6c67a9bd08f97"}, + {file = "SQLAlchemy-1.4.32-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:dd93162615870c976dba43963a24bb418b28448fef584f30755990c134a06a55"}, + {file = "SQLAlchemy-1.4.32-cp38-cp38-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:5a2e73508f939175363d8a4be9dcdc84cf16a92578d7fa86e6e4ca0e6b3667b2"}, + {file = "SQLAlchemy-1.4.32-cp38-cp38-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:bfec934aac7f9fa95fc82147a4ba5db0a8bdc4ebf1e33b585ab8860beb10232f"}, + {file = "SQLAlchemy-1.4.32-cp38-cp38-win32.whl", hash = "sha256:bb42f9b259c33662c6a9b866012f6908a91731a419e69304e1261ba3ab87b8d1"}, + {file = "SQLAlchemy-1.4.32-cp38-cp38-win_amd64.whl", hash = "sha256:7ff72b3cc9242d1a1c9b84bd945907bf174d74fc2519efe6184d6390a8df478b"}, + {file = "SQLAlchemy-1.4.32-cp39-cp39-macosx_10_15_x86_64.whl", hash = "sha256:5dc9801ae9884e822ba942ca493642fb50f049c06b6dbe3178691fce48ceb089"}, + {file = "SQLAlchemy-1.4.32-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:e4607d2d16330757818c9d6fba322c2e80b4b112ff24295d1343a80b876eb0ed"}, + {file = "SQLAlchemy-1.4.32-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:20e9eba7fd86ef52e0df25bea83b8b518dfdf0bce09b336cfe51671f52aaaa3f"}, + {file = "SQLAlchemy-1.4.32-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:290cbdf19129ae520d4bdce392648c6fcdbee763bc8f750b53a5ab51880cb9c9"}, + {file = "SQLAlchemy-1.4.32-cp39-cp39-win32.whl", hash = "sha256:1bbac3e8293b34c4403d297e21e8f10d2a57756b75cff101dc62186adec725f5"}, + {file = "SQLAlchemy-1.4.32-cp39-cp39-win_amd64.whl", hash = "sha256:b3f1d9b3aa09ab9adc7f8c4b40fc3e081eb903054c9a6f9ae1633fe15ae503b4"}, + {file = "SQLAlchemy-1.4.32.tar.gz", hash = "sha256:6fdd2dc5931daab778c2b65b03df6ae68376e028a3098eb624d0909d999885bc"}, ] stack-data = [ {file = "stack_data-0.2.0-py3-none-any.whl", hash = "sha256:999762f9c3132308789affa03e9271bbbe947bf78311851f4d485d8402ed858e"}, @@ -2207,8 +2135,8 @@ urllib3 = [ {file = "urllib3-1.26.8.tar.gz", hash = "sha256:0e7c33d9a63e7ddfcb86780aac87befc2fbddf46c58dbb487e0855f7ceec283c"}, ] virtualenv = [ - {file = "virtualenv-20.13.2-py2.py3-none-any.whl", hash = "sha256:e7b34c9474e6476ee208c43a4d9ac1510b041c68347eabfe9a9ea0c86aa0a46b"}, - {file = "virtualenv-20.13.2.tar.gz", hash = "sha256:01f5f80744d24a3743ce61858123488e91cb2dd1d3bdf92adaf1bba39ffdedf0"}, + {file = "virtualenv-20.13.3-py2.py3-none-any.whl", hash = "sha256:dd448d1ded9f14d1a4bfa6bfc0c5b96ae3be3f2d6c6c159b23ddcfd701baa021"}, + {file = "virtualenv-20.13.3.tar.gz", hash = "sha256:e9dd1a1359d70137559034c0f5433b34caf504af2dc756367be86a5a32967134"}, ] wcwidth = [ {file = "wcwidth-0.2.5-py2.py3-none-any.whl", hash = "sha256:beb4802a9cebb9144e99086eff703a642a13d6a0052920003a230f3294bbe784"}, diff --git a/pyproject.toml b/pyproject.toml index 8d3dda3f2..1a2f29375 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -76,13 +76,15 @@ ekonav = "ekomark.navigator:launch_navigator" [tool.poe.tasks] test = "pytest tests" bench = "pytest benchmarks" -bench-iso = "pytest benchmarks -m isolated" -bench-run = "pytest benchmarks -m 'not isolated'" +bench-iso.cmd = "pytest benchmarks -m isolated" +bench-iso.env.NUMBA_DISABLE_JIT.default = "0" +bench-run.cmd = "pytest benchmarks -m 'not isolated'" +bench-run.env.NUMBA_DISABLE_JIT.default = "0" lint = "pylint src/**/*.py -E" lint-warnings = "pylint src/**/*.py --exit-zero" sandbox = "python benchmarks/sandbox.py" -nav = "ekonav" -navigator = "ekonav" +nav = "ekonav --config benchmarks/banana.yaml" +navigator = "ekonav --config benchmarks/banana.yaml" docs = { "shell" = "cd doc; make html" } docs-view = { "shell" = "cd doc; make view" } docs-server = { "shell" = "cd doc; make server" } @@ -101,7 +103,7 @@ addopts = [ '--cov-report=xml', '--strict-markers', ] -env = ["NUMBA_DISABLE_JIT=1"] +env = ["D:NUMBA_DISABLE_JIT=1"] markers = ["isolated: marks benchmarks as isolated"] [tool.pylint.master] diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index f42b972b6..ac18eb7b4 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -6,14 +6,17 @@ """ import logging +import os import time +from multiprocessing import Pool import numba as nb import numpy as np from scipy import integrate from .. import anomalous_dimensions as ad -from .. import beta, interpolation, mellin +from .. import interpolation, mellin +from .. import scale_variations as sv from ..basis_rotation import full_labels, singlet_labels from ..kernels import non_singlet as ns from ..kernels import singlet as s @@ -21,83 +24,103 @@ logger = logging.getLogger(__name__) +sv_mode_dict = dict( + zip( + [None, "exponentiated", "expanded"], + [sv.unvaried, sv.mode_exponentiated, sv.mode_expanded], + ) +) + -@nb.njit("c16[:](u1,string,c16,u1,f8)", cache=True) -def gamma_ns_fact(order, mode, n, nf, L): +@nb.njit("c16(c16[:,:],string)") +def select_singlet_element(ker, mode): """ - Adjust the anomalous dimensions with the scale variations. + Select element of the singlet matrix Parameters ---------- - order : int - perturbation order mode : str sector element - n : complex - Melling moment - nf : int - number of active flavors - L : float - logarithmic ratio of factorization and renormalization scale + ker : numpy.ndarray + singlet integration kernel Returns ------- - gamma_ns : numpy.ndarray - adjusted non-singlet anomalous dimensions + ker : complex + singlet integration kernel element """ - gamma_ns = ad.gamma_ns(order, mode[-1], n, nf) - # since we are modifying *in-place* be carefull, that the order matters! - # and indeed, we need to adjust the high elements first - if order >= 2: - gamma_ns[2] -= ( - 2 * beta.beta(0, nf) * gamma_ns[1] * L - + (beta.beta(1, nf) * L - beta.beta(0, nf) ** 2 * L**2) * gamma_ns[0] - ) - if order >= 1: - gamma_ns[1] -= beta.beta(0, nf) * gamma_ns[0] * L - return gamma_ns + k = 0 if mode[2] == "q" else 1 + l = 0 if mode[3] == "q" else 1 + return ker[k, l] + + +spec = [ + ("is_singlet", nb.boolean), + ("is_log", nb.boolean), + ("logx", nb.float64), + ("u", nb.float64), +] -@nb.njit("c16[:,:,:](u1,c16,u1,f8)", cache=True) -def gamma_singlet_fact(order, n, nf, L): +@nb.experimental.jitclass(spec) +class QuadKerBase: """ - Adjust the anomalous dimensions with the scale variations. + Manage the common part of Mellin inversion integral Parameters ---------- - order : int - perturbation order + u : float + quad argument + is_log : boolean + is a logarithmic interpolation + logx : float + Mellin inversion point mode : str sector element - n : complex - Melling moment - nf : int - number of active flavors - L : float - logarithmic ratio of factorization and renormalization scale - - Returns - ------- - gamma_singlet : numpy.ndarray - adjusted singlet anomalous dimensions """ - gamma_singlet = ad.gamma_singlet( - order, - n, - nf, - ) - # concerning order: see comment at gamma_ns_fact - if order >= 2: - gamma_singlet[2] -= ( - 2 * beta.beta(0, nf) * gamma_singlet[1] * L - + (beta.beta(1, nf) * L - beta.beta(0, nf) ** 2 * L**2) * gamma_singlet[0] - ) - if order >= 1: - gamma_singlet[1] -= beta.beta(0, nf) * gamma_singlet[0] * L - return gamma_singlet + def __init__(self, u, is_log, logx, mode): + self.is_singlet = mode[0] == "S" + self.is_log = is_log + self.u = u + self.logx = logx + + @property + def path(self): + """Returns the associated instance of :class:`eko.mellin.Path`""" + return mellin.Path(self.u, self.logx, self.is_singlet) + + @property + def n(self): + """Returns the Mellin moment N""" + return self.path.n + + def integrand( + self, + areas, + ): + """ + Get transformation to Mellin space integral -@nb.njit("f8(f8,u1,string,string,b1,f8,f8[:,:],f8,f8,f8,f8,u4,u1)", cache=True) + Parameters + ---------- + areas : tuple + basis function configuration + + Returns + ------- + base_integrand: complex + common mellin inversion intgrand + """ + if self.logx == 0.0: + return 0.0 + pj = interpolation.evaluate_grid(self.path.n, self.is_log, self.logx, areas) + if pj == 0.0: + return 0.0 + return self.path.prefactor * pj * self.path.jac + + +@nb.njit("f8(f8,u1,string,string,b1,f8,f8[:,:],f8,f8,f8,f8,u4,u1,u1)", cache=True) def quad_ker( u, order, @@ -112,9 +135,10 @@ def quad_ker( L, ev_op_iterations, ev_op_max_order, + sv_mode, ): """ - Raw kernel inside quad. + Raw evolution kernel inside quad. Parameters ---------- @@ -144,42 +168,40 @@ def quad_ker( number of evolution steps ev_op_max_order : int perturbative expansion order of U + sv_mode: int + use scale variation mode 0: none, 1: exponentiated, 2: expanded Returns ------- ker : float evaluated integration kernel """ - is_singlet = mode[0] == "S" - # get transformation to N integral - if logx == 0.0: - return 0.0 - r = 0.4 * 16.0 / (-logx) - if is_singlet: - o = 1.0 - else: - o = 0.0 - n = mellin.Talbot_path(u, r, o) - jac = mellin.Talbot_jac(u, r, o) - # check PDF is active - if is_log: - pj = interpolation.log_evaluate_Nx(n, logx, areas) - else: - pj = interpolation.evaluate_Nx(n, logx, areas) - if pj == 0.0: + ker_base = QuadKerBase(u, is_log, logx, mode) + integrand = ker_base.integrand(areas) + if integrand == 0.0: return 0.0 + # compute the actual evolution kernel - if is_singlet: - gamma_singlet = gamma_singlet_fact(order, n, nf, L) + if ker_base.is_singlet: + gamma_singlet = ad.gamma_singlet(order, ker_base.n, nf) + # scale var A is directly applied on gamma + if sv_mode == sv.mode_exponentiated: + gamma_singlet = sv.exponentiated.gamma_variation( + gamma_singlet, order, nf, L + ) ker = s.dispatcher( order, method, gamma_singlet, a1, a0, nf, ev_op_iterations, ev_op_max_order ) - # select element of matrix - k = 0 if mode[2] == "q" else 1 - l = 0 if mode[3] == "q" else 1 - ker = ker[k, l] + # scale var B is applied on the kernel + if sv_mode == sv.mode_expanded: + ker = np.ascontiguousarray(ker) @ np.ascontiguousarray( + sv.expanded.singlet_variation(gamma_singlet, a1, order, nf, L) + ) + ker = select_singlet_element(ker, mode) else: - gamma_ns = gamma_ns_fact(order, mode, n, nf, L) + gamma_ns = ad.gamma_ns(order, mode[-1], ker_base.n, nf) + if sv_mode == sv.mode_exponentiated: + gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) ker = ns.dispatcher( order, method, @@ -189,9 +211,11 @@ def quad_ker( nf, ev_op_iterations, ) + if sv_mode == sv.mode_expanded: + ker = ker * sv.expanded.non_singlet_variation(gamma_ns, a1, order, nf, L) + # recombine everthing - mellin_prefactor = complex(0.0, -1.0 / np.pi) - return np.real(mellin_prefactor * ker * pj * jac) + return np.real(ker * integrand) class Operator: @@ -226,6 +250,37 @@ def __init__(self, config, managers, nf, q2_from, q2_to, mellin_cut=5e-2): self._mellin_cut = mellin_cut self.op_members = {} + @property + def fact_to_ren(self): + r"""Returns the factor :math:`(\mu_F/\mu_R)^2`""" + return self.config["fact_to_ren"] + + @property + def sv_mode(self): + """Returns the scale variation mode""" + return sv_mode_dict[self.config["ModSV"]] + + @property + def int_disp(self): + """Returns the interpolation dispatcher""" + return self.managers["interpol_dispatcher"] + + @property + def grid_size(self): + """Returns the grid size""" + return self.int_disp.xgrid.size + + @property + def a_s(self): + """Returns the computed values for :math:`a_s`""" + sc = self.managers["strong_coupling"] + a0 = sc.a_s( + self.q2_from / self.fact_to_ren, fact_scale=self.q2_from, nf_to=self.nf + ) + a1 = sc.a_s(self.q2_to / self.fact_to_ren, fact_scale=self.q2_to, nf_to=self.nf) + return (a0, a1) + + @property def labels(self): """ Compute necessary sector labels to compute. @@ -254,18 +309,19 @@ def labels(self): labels.extend(singlet_labels) return labels - def compute(self): - """compute the actual operators (i.e. run the integrations)""" - # Generic parameters - int_disp = self.managers["interpol_dispatcher"] - grid_size = len(int_disp.xgrid) - - # init all ops with identity or zeros if we skip them - labels = self.labels() - eye = OpMember(np.eye(grid_size), np.zeros((grid_size, grid_size))) - zero = OpMember(*[np.zeros((grid_size, grid_size))] * 2) + @property + def quad_ker(self): + """Integrand function""" + return quad_ker + + def initialize_op_members(self): + """Init all ops with identity or zeros if we skip them""" + eye = OpMember( + np.eye(self.grid_size), np.zeros((self.grid_size, self.grid_size)) + ) + zero = OpMember(*[np.zeros((self.grid_size, self.grid_size))] * 2) for n in full_labels: - if n in labels: + if n in self.labels: # off diag singlet are zero if n in ["S_qg", "S_gq"]: self.op_members[n] = zero.copy() @@ -273,17 +329,78 @@ def compute(self): self.op_members[n] = eye.copy() else: self.op_members[n] = zero.copy() - # skip computation + + def run_op_integration( + self, + log_grid, + ): + """ + Run the integration for each grid point + + Parameters + ---------- + log_grid : tuple(k, logx) + log grid point with relative index + + Returns + ------- + column : list + computed operators at the give grid point + + """ + column = [] + k, logx = log_grid + start_time = time.perf_counter() + # iterate basis functions + for l, bf in enumerate(self.int_disp): + if k == l and l == self.grid_size - 1: + continue + temp_dict = {} + # iterate sectors + for label in self.labels: + res = integrate.quad( + self.quad_ker, + 0.5, + 1.0 - self._mellin_cut, + args=( + self.config["order"], + label, + self.config["method"], + self.int_disp.log, + logx, + bf.areas_representation, + self.a_s[1], + self.a_s[0], + self.nf, + np.log(self.fact_to_ren), + self.config["ev_op_iterations"], + self.config["ev_op_max_order"], + self.sv_mode, + ), + epsabs=1e-12, + epsrel=1e-5, + limit=100, + full_output=1, + ) + temp_dict[label] = res[:2] + column.append(temp_dict) + + print( + f"Evolution: computing operators: - {k+1}/{self.grid_size} took: {(time.perf_counter() - start_time):6f} s" # pylint: disable=line-too-long + ) + return column + + def compute(self): + """compute the actual operators (i.e. run the integrations)""" + self.initialize_op_members() + + # skip computation ? if np.isclose(self.q2_from, self.q2_to): logger.info("Evolution: skipping unity operator at %e", self.q2_from) self.copy_ns_ops() return + tot_start_time = time.perf_counter() - # setup ingredients - sc = self.managers["strong_coupling"] - fact_to_ren = self.config["fact_to_ren"] - a0 = sc.a_s(self.q2_from / fact_to_ren, fact_scale=self.q2_from, nf_to=self.nf) - a1 = sc.a_s(self.q2_to / fact_to_ren, fact_scale=self.q2_to, nf_to=self.nf) logger.info( "Evolution: computing operators %e -> %e, nf=%d", self.q2_from, @@ -292,60 +409,36 @@ def compute(self): ) logger.info( "Evolution: µ_R^2 distance: %e -> %e", - self.q2_from / fact_to_ren, - self.q2_to / fact_to_ren, + self.q2_from / self.fact_to_ren, + self.q2_to / self.fact_to_ren, ) - logger.info("Evolution: a_s distance: %e -> %e", a0, a1) + if self.sv_mode != 0: + logger.info( + "Scale Variation: (µ_F/µ_R)^2 = %e, mode: %s", + self.fact_to_ren, + "exponentiated" if self.sv_mode == 1 else "expanded", + ) + logger.info("Evolution: a_s distance: %e -> %e", self.a_s[0], self.a_s[1]) logger.info( "Evolution: order: %d, solution strategy: %s", self.config["order"], self.config["method"], ) - logger.info("Evolution: computing operators - 0/%d", grid_size) - # iterate output grid - for k, logx in enumerate(np.log(int_disp.xgrid_raw)): - start_time = time.perf_counter() - # iterate basis functions - for l, bf in enumerate(int_disp): - if k == l and l == grid_size - 1: - continue - # iterate sectors - for label in labels: - # compute and set - res = integrate.quad( - quad_ker, - 0.5, - 1.0 - self._mellin_cut, - args=( - self.config["order"], - label, - self.config["method"], - int_disp.log, - logx, - bf.areas_representation, - a1, - a0, - self.nf, - np.log(fact_to_ren), - self.config["ev_op_iterations"], - self.config["ev_op_max_order"], - ), - epsabs=1e-12, - epsrel=1e-5, - limit=100, - full_output=1, - ) - val, err = res[:2] - self.op_members[label].value[k][l] = val - self.op_members[label].error[k][l] = err - logger.info( - "Evolution: computing operators - %d/%d took: %f s", - k + 1, - grid_size, - time.perf_counter() - start_time, + # run integration in parallel for each grid point + with Pool(int(os.cpu_count() / 2)) as pool: + res = pool.map( + self.run_op_integration, + enumerate(np.log(self.int_disp.xgrid_raw)), ) + # collect results + for k, row in enumerate(res): + for l, entry in enumerate(row): + for label, (val, err) in entry.items(): + self.op_members[label].value[k][l] = val + self.op_members[label].error[k][l] = err + # closing comment logger.info("Evolution: Total time %f s", time.perf_counter() - tot_start_time) # copy non-singlet kernels, if necessary diff --git a/src/eko/evolution_operator/grid.py b/src/eko/evolution_operator/grid.py index aa9e80b76..bea651f1c 100644 --- a/src/eko/evolution_operator/grid.py +++ b/src/eko/evolution_operator/grid.py @@ -11,6 +11,8 @@ import numpy as np +from eko.thresholds import flavor_shift, is_downward_path + from .. import basis_rotation as br from .. import matching_conditions, member from ..matching_conditions.operator_matrix_element import OperatorMatrixElement @@ -125,6 +127,7 @@ def from_dict( config["debug_skip_singlet"] = operators_card["debug_skip_singlet"] config["debug_skip_non_singlet"] = operators_card["debug_skip_non_singlet"] config["HQ"] = theory_card["HQ"] + config["ModSV"] = theory_card["ModSV"] q2_grid = np.array(operators_card["Q2grid"], np.float_) intrinsic_range = [] if int(theory_card["IC"]) == 1: @@ -153,7 +156,7 @@ def get_threshold_operators(self, path): It computes and stores the necessary macthing operators Parameters ---------- - path: list(PathSegment) + path: list(`eko.thresholds.PathSegment`) thresholds path Returns @@ -162,9 +165,12 @@ def get_threshold_operators(self, path): """ # The base area is always that of the reference q thr_ops = [] + # is_downward point to smaller nf + is_downward = is_downward_path(path) + shift = flavor_shift(is_downward) for seg in path[:-1]: new_op_key = seg.tuple - ome = OperatorMatrixElement(self.config, self.managers, seg.is_backward) + ome = OperatorMatrixElement(self.config, self.managers, is_downward) if new_op_key not in self._threshold_operators: # Compute the operator and store it logger.info("Prepare threshold operator") @@ -178,10 +184,13 @@ def get_threshold_operators(self, path): # Compute the matching conditions and store it if seg.q2_to not in self._matching_operators: thr_config = self.managers["thresholds_config"] - # is_backward point to the smaller q2 - shift = 3 if not seg.is_backward else 4 kthr = thr_config.thresholds_ratios[seg.nf - shift] - ome.compute(seg.q2_to, np.log(kthr), self.config["HQ"] == "MSBAR") + ome.compute( + seg.q2_to, + seg.nf - shift + 3, + np.log(kthr), + self.config["HQ"] == "MSBAR", + ) self._matching_operators[seg.q2_to] = ome.ome_members return thr_ops @@ -213,7 +222,7 @@ def compute(self, q2grid=None): def generate(self, q2): """ - Computes an single EKO. + Computes a single EKO. Parameters ---------- @@ -235,7 +244,8 @@ def generate(self, q2): ) operator.compute() intrinsic_range = self.config["intrinsic_range"] - if path[-1].is_backward: + is_downward = is_downward_path(path) + if is_downward: intrinsic_range = [4, 5, 6] final_op = physical.PhysicalOperator.ad_to_evol_map( operator.op_members, @@ -243,16 +253,14 @@ def generate(self, q2): operator.q2_to, intrinsic_range, ) - # and multiply the lower ones from the right - for i, op in reversed(list(enumerate(thr_ops))): - is_backward = path[i].is_backward + 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 ) # join with the basis rotation, since matching requires c+ (or likewise) - if is_backward: + if is_downward: matching = matching_conditions.MatchingCondition.split_ad_to_evol_map( self._matching_operators[op.q2_to], op.nf - 1, @@ -281,7 +289,7 @@ def generate(self, q2): "operators": values, "operator_errors": errors, "alphas": self.managers["strong_coupling"].a_s( - q2 / fact_to_ren, fact_scale=q2 + q2 / fact_to_ren, fact_scale=q2, nf_to=path[-1].nf ) * 4.0 * np.pi, diff --git a/src/eko/interpolation.py b/src/eko/interpolation.py index 442bbd847..bb0ade7b0 100644 --- a/src/eko/interpolation.py +++ b/src/eko/interpolation.py @@ -176,6 +176,34 @@ def evaluate_Nx(N, logx, area_list): return res +@nb.njit("c16(c16,b1,f8,f8[:,:])", cache=True) +def evaluate_grid(N, is_log, logx, area_list): + """ + Evaluate interpolator in N-space + + Parameters + ---------- + N : complex + Mellin variable + is_log : boolean + is a logarithmic interpolation + logx : float + logarithm of inversion point + area_list : list + area configuration of basis function + + Returns + ------- + pj : float + basis function * inversion factor + """ "" + if is_log: + pj = log_evaluate_Nx(N, logx, area_list) + else: + pj = evaluate_Nx(N, logx, area_list) + return pj + + # TODO lift to runcard? _atol_eps = 10 * np.finfo(float).eps @@ -482,19 +510,15 @@ def from_dict(cls, operators_card, mode_N=True): * - Name - Type - - default - description * - ``interpolation_xgrid`` - :py:obj:`list(float)` - - [required] - the interpolation grid * - ``interpolation_polynomial_degree`` - :py:obj:`int` - - ``4`` - polynomial degree of the interpolating function * - ``interpolation_is_log`` - :py:obj:`bool` - - ``True`` - use logarithmic interpolation? Parameters diff --git a/src/eko/kernels/singlet.py b/src/eko/kernels/singlet.py index a55fac2c8..aa5d19f33 100644 --- a/src/eko/kernels/singlet.py +++ b/src/eko/kernels/singlet.py @@ -532,12 +532,9 @@ def dispatcher( # pylint: disable=too-many-return-statements elif method == "decompose-exact": if order == 1: return nlo_decompose_exact(gamma_singlet, a1, a0, nf) - else: - return nnlo_decompose_exact(gamma_singlet, a1, a0, nf) + return nnlo_decompose_exact(gamma_singlet, a1, a0, nf) elif method == "decompose-expanded": if order == 1: return nlo_decompose_expanded(gamma_singlet, a1, a0, nf) - else: - return nnlo_decompose_expanded(gamma_singlet, a1, a0, nf) - else: - raise NotImplementedError("Selected method is not implemented") + return nnlo_decompose_expanded(gamma_singlet, a1, a0, nf) + raise NotImplementedError("Selected method is not implemented") diff --git a/src/eko/matching_conditions/__init__.py b/src/eko/matching_conditions/__init__.py index 3b7ae3c72..96cf0356e 100644 --- a/src/eko/matching_conditions/__init__.py +++ b/src/eko/matching_conditions/__init__.py @@ -39,8 +39,6 @@ def split_ad_to_evol_map( threshold value intrinsic_range : list list of intrinsic quark pids - is_backward: bool - True for backward evolution """ m = { diff --git a/src/eko/matching_conditions/operator_matrix_element.py b/src/eko/matching_conditions/operator_matrix_element.py index 708ca218e..b610f10a0 100644 --- a/src/eko/matching_conditions/operator_matrix_element.py +++ b/src/eko/matching_conditions/operator_matrix_element.py @@ -286,7 +286,7 @@ def labels(self): # labels.extend(["S_qH"]) return labels - def compute(self, q2, L, is_msbar): + def compute(self, q2, nf, L, is_msbar): """ compute the actual operators (i.e. run the integrations) @@ -294,6 +294,8 @@ def compute(self, q2, L, is_msbar): ---------- q2: float threshold scale + nf: int + number of active flavor below threshold L: float log of K threshold squared is_msbar: bool @@ -319,7 +321,8 @@ def compute(self, q2, L, is_msbar): self.copy_ome() return - a_s = self.sc.a_s(q2 / self.config["fact_to_ren"], q2) + # Note that here you need to use a_s^{nf+1}(q2) + a_s = self.sc.a_s(q2 / self.config["fact_to_ren"], q2, nf + 1) tot_start_time = time.perf_counter() logger.info("Matching: computing operators - 0/%d", grid_size) diff --git a/src/eko/mellin.py b/src/eko/mellin.py index 3abb6fc5d..c7a717590 100644 --- a/src/eko/mellin.py +++ b/src/eko/mellin.py @@ -166,8 +166,7 @@ def edge_path(t, m, c, phi): """ if t < 0.5: # turning point: path is not differentiable in this point return c + (0.5 - t) * m * np.exp(complex(0, -phi)) - else: - return c + (t - 0.5) * m * np.exp(complex(0, +phi)) + return c + (t - 0.5) * m * np.exp(complex(0, +phi)) @nb.njit("c16(f8,f8,f8,f8)", cache=True) @@ -196,5 +195,52 @@ def edge_jac(t, m, _c, phi): """ if t < 0.5: # turning point: jacobian is not continuous here return -m * np.exp(complex(0, -phi)) - else: - return +m * np.exp(complex(0, phi)) + return +m * np.exp(complex(0, phi)) + + +spec = [ + ("t", nb.float64), + ("r", nb.float64), + ("o", nb.int8), +] + + +@nb.experimental.jitclass(spec) +class Path: + """ + Mellin path dispatcher + + Parameters + ---------- + t : float + way parameter + logx : float + Mellin inversion point + axis_offset: bool + add offset on the real axis + """ + + def __init__(self, t, logx, axis_offset): + self.t = t + # TODO: shall we use: 0.4 * 16.0 / (1.0 - logx) ? + self.r = 0.4 * 16.0 / (-logx) + if axis_offset: + self.o = 1.0 + else: + self.o = 0.0 + + # TODO: make also the other 2 paths available ?? + @property + def n(self): + """Returns the Mellin moment N""" + return Talbot_path(self.t, self.r, self.o) + + @property + def jac(self): + """Returns the Jacobian of the Mellin path""" + return Talbot_jac(self.t, self.r, self.o) + + @property + def prefactor(self): + r"""Returns the mellin inversion prefactor :math:`-\frac{i}{\pi}`""" + return complex(0.0, -1.0 / np.pi) diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index 6c3813d63..9b984704d 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -10,6 +10,7 @@ from .beta import b, beta from .gamma import gamma from .strong_coupling import StrongCoupling, invert_matching_coeffs +from .thresholds import ThresholdsAtlas, flavor_shift, is_downward_path def ker_exact(a0, a1, order, nf): @@ -163,12 +164,12 @@ def ker_dispatcher(q2_to, q2m_ref, strong_coupling, fact_to_ren, nf): ker: Expanded or exact |MSbar| kernel """ - a0 = strong_coupling.a_s(q2m_ref / fact_to_ren, q2m_ref) - a1 = strong_coupling.a_s(q2_to / fact_to_ren, q2_to) + a0 = strong_coupling.a_s(q2m_ref / fact_to_ren, q2m_ref, nf) + a1 = strong_coupling.a_s(q2_to / fact_to_ren, q2_to, nf) method = strong_coupling.method order = strong_coupling.order if method == "expanded": - return ker_expanded(a0, a1, order, nf) + return ker_expanded(a0, float(a1), order, nf) return ker_exact(a0, a1, order, nf) @@ -267,11 +268,7 @@ def rge(m2, q2m_ref, strong_coupling, fact_to_ren, nf_ref): def evolve( - m2_ref, - q2m_ref, - strong_coupling, - fact_to_ren, - q2_to, + m2_ref, q2m_ref, strong_coupling, fact_to_ren, q2_to, nf_ref=None, nf_to=None ): r""" Perform the |MSbar| mass evolution up to given scale. @@ -288,55 +285,57 @@ def evolve( any q fact_to_ren: float :math:`\mu_F^2/\mu_R^2` - q2_to: float, optional + q2_to: float scale at which the mass is computed + nf_ref: int + number of flavor active at the reference scale + nf_to: int + number of flavor active at the target scale Returns ------- m2 : float :math:`m_{\overline{MS}}(\mu_2)^2` """ - # evolution might involve different number of active flavors. - # Find out the evolution path (always sorted) - q2_low, q2_high = sorted([q2m_ref, q2_to]) - - area_walls = np.array(strong_coupling.thresholds.area_walls) - path = np.concatenate( - ( - [q2_low], - area_walls[(q2_low < area_walls) & (area_walls < q2_high)], - [q2_high], - ) + thr_atlas = ThresholdsAtlas( + np.array(strong_coupling.thresholds.area_walls)[1:-1], + q2m_ref, + nf_ref, + strong_coupling.thresholds.thresholds_ratios, ) - nf_init = len(area_walls[q2_low >= area_walls]) + 2 - nf_final = len(area_walls[q2_high > area_walls]) + 2 + path = thr_atlas.path(q2_to, nf_to) + is_downward = is_downward_path(path) + shift = flavor_shift(is_downward) ev_mass = 1.0 - is_downward_path = bool(q2m_ref > q2_to) - for i, nf in enumerate(np.arange(nf_init, nf_final + 1)): - q2_init = path[i] - q2_final = path[i + 1] - # if you are going backward - # need to reverse the evolution in each path segment - if is_downward_path: - m_coeffs = compute_matching_coeffs_down(nf - 1) - q2_init, q2_final = q2_final, q2_init - shift = 4 - else: - m_coeffs = compute_matching_coeffs_up(nf) - shift = 3 - fact = 1.0 - # shift - for pto in range(1, strong_coupling.order + 1): - for l in range(pto + 1): - as_thr = strong_coupling.a_s(q2_final / fact_to_ren, q2_final) - # TODO: do we need to add np.log(fac_to_ren) here ??? - L = np.log(strong_coupling.thresholds.thresholds_ratios[pto - shift]) - fact += as_thr**pto * L**l * m_coeffs[pto, l] - ev_mass *= ( - fact - * ker_dispatcher(q2_final, q2_init, strong_coupling, fact_to_ren, nf) ** 2 - ) + for k, seg in enumerate(path): + # skip a very short segment, but keep the matching + ker_evol = 1.0 + if not np.isclose(seg.q2_from, seg.q2_to): + ker_evol = ( + ker_dispatcher( + seg.q2_to, seg.q2_from, strong_coupling, fact_to_ren, seg.nf + ) + ** 2 + ) + # apply matching condition + if k < len(path) - 1: + # TODO: do we need to add np.log(fac_to_ren) here ??? + L = np.log(thr_atlas.thresholds_ratios[seg.nf - shift]) + m_coeffs = ( + compute_matching_coeffs_down(seg.nf - 1) + if is_downward + else compute_matching_coeffs_up(seg.nf) + ) + matching = 1.0 + for pto in range(1, strong_coupling.order + 1): + for l in range(pto + 1): + as_thr = strong_coupling.a_s( + seg.q2_to / fact_to_ren, seg.q2_to, seg.nf - shift + 4 + ) + matching += as_thr**pto * L**l * m_coeffs[pto, l] + ker_evol *= matching + ev_mass *= ker_evol return m2_ref * ev_mass @@ -421,7 +420,8 @@ def sc(thr_masses): # you need to evolve it until nf_target patch wall is reached: # for backward you reach the higher, for forward the lower. # len(masses[q2m_ref > masses]) + 3 is the nf at the given reference scale - if nf_target != len(masses[q2m_ref > masses]) + 3: + nf_ref = len(masses[q2m_ref > masses]) + 3 + if nf_target != nf_ref: q2_to = masses[q_idx + shift] m2_ref = evolve( m2_ref, @@ -429,6 +429,8 @@ def sc(thr_masses): sc(masses), fact_to_ren, q2_to, + nf_ref=nf_ref, + nf_to=nf_target, ) q2m_ref = q2_to diff --git a/src/eko/scale_variations/__init__.py b/src/eko/scale_variations/__init__.py new file mode 100644 index 000000000..942988e77 --- /dev/null +++ b/src/eko/scale_variations/__init__.py @@ -0,0 +1,11 @@ +# -*- coding: utf-8 -*- +""" +This module contain the possible scale variations integrals. +""" + +from . import expanded, exponentiated + +unvaried = 0 +mode_exponentiated = 1 +mode_expanded = 2 +"""Scale Variation modes""" diff --git a/src/eko/scale_variations/expanded.py b/src/eko/scale_variations/expanded.py new file mode 100644 index 000000000..610f2f225 --- /dev/null +++ b/src/eko/scale_variations/expanded.py @@ -0,0 +1,173 @@ +# -*- coding: utf-8 -*- +r""" +This module contains the scale variation operator in ``ModSV=expanded`` +""" + + +import numba as nb +import numpy as np + +from .. import beta + + +@nb.njit(["c16(c16[:],f8)", "c16[:,:](c16[:,:,:],f8)"], cache=True) +def gamma_1_variation(gamma, L): + r""" + Computes the |NLO| anomalous dimension variation. + + Parameters + ---------- + gamma : numpy.ndarray + anomalous dimensions + L : float + logarithmic ratio of factorization and renormalization scale + + Returns + ------- + gamma_1 : complex + variation to :math:`\gamma^{(1)}` + """ + return -L * gamma[0] + + +@nb.njit(["c16(c16[:],f8,f8,c16)", "c16[:,:](c16[:,:,:],f8,f8,c16[:,:])"], cache=True) +def gamma_2_variation(gamma, L, beta0, g0e2): + r""" + Computes the |NNLO| anomalous dimension variation. + + Parameters + ---------- + gamma : numpy.ndarray + anomalous dimensions + L : float + logarithmic ratio of factorization and renormalization scale + beta0: float + :math:`\beta_0` + g0e2: complex + :math:`\gamma^{(0),2}` + + Returns + ------- + gamma_2 : complex + variation to :math:`\gamma^{(2)}` + """ + return -gamma[1] * L + 1 / 2 * (beta0 * gamma[0] + g0e2) * L**2 + + +@nb.njit( + [ + "c16(c16[:],f8,f8,f8,c16,c16,c16)", + "c16[:,:](c16[:,:,:],f8,f8,f8,c16[:,:],c16[:,:],c16[:,:])", + ], + cache=True, +) +def gamma_3_variation(gamma, L, beta0, beta1, g0e2, g0e3, g1g0): + r""" + Computes the |N3LO| anomalous dimension variation. + + Parameters + ---------- + gamma : numpy.ndarray + anomalous dimensions + L : float + logarithmic ratio of factorization and renormalization scale + beta0: float + :math:`\beta_0` + beta0: float + :math:`\beta_1` + g0e2: complex + :math:`\gamma^{(0),2}` + g0e3: complex + :math:`\gamma^{(0),3}` + g1g0: complex + :math:`\gamma^{(1)} \gamma^{(0)}` + + Returns + ------- + gamma_3 : complex + variation to :math:`\gamma^{(3)}` + """ + return ( + -gamma[2] * L + + (1 / 2) * (beta1 * gamma[0] + 2 * beta0 * gamma[1] + 2 * g1g0) * L**2 + - (1 / 6) * (2 * beta0**2 * gamma[0] + 3 * beta0 * g0e2 + g0e3) * L**3 + ) + + +@nb.njit("c16(c16[:],f8,u1,u1,f8)", cache=True) +def non_singlet_variation(gamma, a_s, order, nf, L): + """ + Scale Variation non singlet dispatcher + + Parameters + ---------- + gamma : numpy.ndarray + anomalous dimensions + a_s : float + target coupling value + order : int + perturbation order + nf : int + number of active flavors + L : float + logarithmic ratio of factorization and renormalization scale + + Returns + ------- + sv_ker : numpy.ndarray + scale varion kernel + """ + sv_ker = 1.0 + if order >= 1: + sv_ker += a_s * gamma_1_variation(gamma, L) + if order >= 2: + beta0 = beta.beta_0(nf) + sv_ker += a_s**2 * gamma_2_variation(gamma, L, beta0, gamma[0] ** 2) + if order >= 3: + beta1 = beta.beta(1, nf) + sv_ker += a_s**3 * gamma_3_variation( + gamma, L, beta0, beta1, gamma[0] ** 2, gamma[0] ** 3, gamma[0] * gamma[1] + ) + return sv_ker + + +@nb.njit("c16[:,:](c16[:,:,:],f8,u1,u1,f8)", cache=True) +def singlet_variation(gamma, a_s, order, nf, L): + """ + Scale Variation singlet dispatcher + + Parameters + ---------- + gamma : numpy.ndarray + anomalous dimensions + a_s : float + target coupling value + order : int + perturbation order + nf : int + number of active flavors + L : float + logarithmic ratio of factorization and renormalization scale + + Returns + ------- + sv_ker : numpy.ndarray + scale varion kernel + """ + sv_ker = np.eye(2, dtype=np.complex_) + gamma = np.ascontiguousarray(gamma) + if order >= 1: + sv_ker += a_s * gamma_1_variation(gamma, L) + if order >= 2: + beta0 = beta.beta_0(nf) + gamma0e2 = gamma[0] @ gamma[0] + sv_ker += a_s**2 * gamma_2_variation(gamma, L, beta0, gamma0e2) + if order >= 3: + beta1 = beta.beta(1, nf) + gamma0e3 = gamma0e2 @ gamma[0] + # here the product is not commutative + g1g0 = gamma[1] @ gamma[0] + sv_ker += a_s**3 * gamma_3_variation( + gamma, L, beta0, beta1, gamma0e2, gamma0e3, g1g0 + ) + return sv_ker diff --git a/src/eko/scale_variations/exponentiated.py b/src/eko/scale_variations/exponentiated.py new file mode 100644 index 000000000..d0657de12 --- /dev/null +++ b/src/eko/scale_variations/exponentiated.py @@ -0,0 +1,52 @@ +# -*- coding: utf-8 -*- +r""" +This module contains the scale variation for ``ModSV=exponentiated``. +""" +import numba as nb + +from .. import beta + + +@nb.njit(["c16[:,:,:](c16[:,:,:],u1,u1,f8)", "c16[:](c16[:],u1,u1,f8)"], cache=True) +def gamma_variation(gamma, order, nf, L): + """ + Adjust the anomalous dimensions with the scale variations. + + Parameters + ---------- + gamma : numpy.ndarray + anomalous dimensions + order : int + perturbation order + nf : int + number of active flavors + L : float + logarithmic ratio of factorization and renormalization scale + + Returns + ------- + gamma : numpy.ndarray + adjusted anomalous dimensions + """ + # since we are modifying *in-place* be carefull, that the order matters! + # and indeed, we need to adjust the high elements first + beta0 = beta.beta(0, nf) + beta1 = beta.beta(1, nf) + if order >= 3: + gamma[3] -= ( + 3 * beta0 * L * gamma[2] + + (2 * beta1 * L - 3 * beta0**2 * L**2) * gamma[1] + + ( + beta.beta(2, nf) * L + - 5 / 2 * beta1 * beta0 * L**2 + + beta0**3 * L**3 + ) + * gamma[0] + ) + if order >= 2: + gamma[2] -= ( + 2 * beta0 * gamma[1] * L + (beta1 * L - beta0**2 * L**2) * gamma[0] + ) + if order >= 1: + gamma[1] -= beta0 * gamma[0] * L + return gamma diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index 3e977da7a..0642ecdef 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -151,9 +151,9 @@ class StrongCoupling: scale_ref : float reference scale :math:`\mu_0^2` masses : list(float) - list with quark masses + list with quark masses squared thresholds_ratios : list(float) - list with ratios between the mass and the thresholds + list with ratios between the mass and the thresholds squared order: int Evaluated order of the beta function: ``0`` = LO, ... method : ["expanded", "exact"] @@ -342,7 +342,7 @@ def compute(self, as_ref, nf, scale_from, scale_to): as_new = self.compute_exact(float(as_ref), nf, scale_from, scale_to) else: as_new = as_expanded( - self.order, float(as_ref), nf, scale_from, scale_to + self.order, float(as_ref), nf, scale_from, float(scale_to) ) self.cache[key] = as_new return as_new @@ -371,10 +371,8 @@ def a_s( # Set up the path to follow in order to go from q2_0 to q2_ref final_as = self.as_ref path = self.thresholds.path(scale_to, nf_to) - is_downward_path = False - if len(path) > 1: - is_downward_path = path[1].nf < path[0].nf - shift = 3 if not is_downward_path else 4 + is_downward = thresholds.is_downward_path(path) + shift = thresholds.flavor_shift(is_downward) # as a default assume mu_F^2 = mu_R^2 if fact_scale is None: @@ -394,7 +392,7 @@ def a_s( ) m_coeffs = ( compute_matching_coeffs_down(self.hqm_scheme, seg.nf - 1) - if is_downward_path + if is_downward else compute_matching_coeffs_up(self.hqm_scheme, seg.nf) ) fact = 1.0 @@ -402,7 +400,6 @@ def a_s( for n in range(1, self.order + 1): for l in range(n + 1): fact += new_as**n * L**l * m_coeffs[n, l] - # shift new_as *= fact final_as = new_as return final_as diff --git a/src/eko/thresholds.py b/src/eko/thresholds.py index 36936c7a6..f3f173cb0 100644 --- a/src/eko/thresholds.py +++ b/src/eko/thresholds.py @@ -29,7 +29,7 @@ def __init__(self, q2_from, q2_to, nf): self.nf = nf @property - def is_backward(self): + def is_downward_q2(self): """True if q2_from bigger than q2_to""" return self.q2_from > self.q2_to @@ -260,3 +260,41 @@ def nf(self, q2): """ ref_idx = np.digitize(q2, self.area_walls) return 2 + ref_idx + + +def is_downward_path(path): + """ + Determine if a path is downward: + - in the number of active flavors when the path list contains more than one `PathSegment`, + note this can be different from each `PathSegment.is_downward` + - in :math:`Q^2` when just one single `PathSegment` is given + + Parameters + ---------- + path: list(`eko.thresholds.PathSegment`) + + Returns + ------- + is_downward: bool + True for a downward path + """ + if len(path) == 1: + return path[0].is_downward_q2 + return path[1].nf < path[0].nf + + +def flavor_shift(is_downward): + """ + Determine the shift to number of light flavors + + Parameters + ---------- + is_downward: bool + True for a downward path + + Returns + ------- + shift: 3, 4 + shift to number of light flavors + """ + return 4 if is_downward else 3 diff --git a/src/ekobox/genpdf/cli.py b/src/ekobox/genpdf/cli.py index e60bd259e..02e2c62c8 100644 --- a/src/ekobox/genpdf/cli.py +++ b/src/ekobox/genpdf/cli.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- """ -This module defines the CLI - see :doc:`here ` +This module defines the CLI - see :doc:`here ` """ import click diff --git a/src/ekomark/benchmark/runner.py b/src/ekomark/benchmark/runner.py index aa1f896ef..40ccab8a7 100644 --- a/src/ekomark/benchmark/runner.py +++ b/src/ekomark/benchmark/runner.py @@ -72,7 +72,7 @@ def run_me(self, theory, ocard, _pdf): if self.sandbox: rerun = True ops_id = f"o{ocard['hash'][:6]}_t{theory['hash'][:6]}" - path = f"{banana_cfg.cfg['database_path'].parents[0]}/{ops_id}.yaml" + path = f"{banana_cfg.cfg['paths']['database'].parents[0]}/{ops_id}.tar" if os.path.exists(path): rerun = False @@ -87,8 +87,7 @@ def run_me(self, theory, ocard, _pdf): else: # load print(f"Using cached eko data: {os.path.relpath(path,os.getcwd())}") - with open(path, encoding="utf-8") as o: - out = eko.output.Output.load_yaml(o) + out = eko.output.Output.load_tar(path) if self.plot_operator: diff --git a/src/ekomark/navigator/__init__.py b/src/ekomark/navigator/__init__.py index d9befee26..8197fb35b 100644 --- a/src/ekomark/navigator/__init__.py +++ b/src/ekomark/navigator/__init__.py @@ -2,10 +2,13 @@ """ ekomark specialization of the navigator """ +import argparse +import pathlib + from banana import cfg as banana_cfg from banana import navigator as bnav -from . import navigator +from . import navigator, glob def yelp(*args): @@ -35,21 +38,37 @@ def yelp(*args): return None -h = yelp +def register_globals(configpath): + app = navigator.NavigatorApp(configpath, "sandbox") + glob.app = app -app = navigator.NavigatorApp(banana_cfg.cfg, "sandbox") + glob.glob["yelp"] = yelp + glob.glob["h"] = yelp -# register banana functions -bnav.register_globals(globals(), app) + # register banana functions + bnav.register_globals(glob.glob, glob.app) -# add my functions -dfl = app.log_as_dfd -check_log = app.check_log -plot_pdfs = app.plot_pdfs -display_pdfs = app.display_pdfs -compare = app.compare_external + # add my functions + glob.glob["dfl"] = app.log_as_dfd + glob.glob["check_log"] = app.check_log + glob.glob["plot_pdfs"] = app.plot_pdfs + glob.glob["display_pdfs"] = app.display_pdfs + glob.glob["compare"] = app.compare_external def launch_navigator(): """CLI Entry point""" - return bnav.launch_navigator("eko", "ekomark") + parser = argparse.ArgumentParser() + + parser.add_argument( + "-c", "--config", type=pathlib.Path, default=None, help="Path to config file" + ) + + args = parser.parse_args() + + register_globals(banana_cfg.detect(args.config)) + + # ekomark.navigator makes the globals here (e.g. app, ls, t) available inside IPython + return bnav.launch_navigator( + ["eko", "ekomark", "ekomark.navigator.glob"], skip_cfg=True + ) diff --git a/src/ekomark/navigator/glob.py b/src/ekomark/navigator/glob.py new file mode 100644 index 000000000..893dbb324 --- /dev/null +++ b/src/ekomark/navigator/glob.py @@ -0,0 +1,2 @@ +app = None +glob = globals() diff --git a/tests/eko/test_ev_op_grid.py b/tests/eko/test_ev_op_grid.py index 37fe0ceb0..4354346a7 100644 --- a/tests/eko/test_ev_op_grid.py +++ b/tests/eko/test_ev_op_grid.py @@ -36,6 +36,7 @@ def _get_setup(self, use_FFNS): "MaxNfPdf": 6, "MaxNfAs": 6, "HQ": "POLE", + "ModSV": None, } operators_card = { "Q2grid": [1, 100**2], @@ -134,3 +135,22 @@ def test_alphas(self): ) sv_opg = sv_opgrid.compute(3) assert opg[3]["alphas"] < sv_opg[3]["alphas"] + + def test_mod_expanded(self): + theory_update = { + "PTO": 1, + "ModSV": "expanded", + } + epsilon = 1e-1 + for ffns, nf0 in zip([False, True], [5, 3]): + theory_update["nf0"] = nf0 + opgrid = self._get_operator_grid(use_FFNS=ffns, theory_update=theory_update) + opg = opgrid.compute(3) + theory_update["fact_to_ren_scale_ratio"] = 1.0 + epsilon + sv_opgrid = self._get_operator_grid( + use_FFNS=ffns, theory_update=theory_update + ) + sv_opg = sv_opgrid.compute(3) + np.testing.assert_allclose( + opg[3]["operators"], sv_opg[3]["operators"], atol=0.7 * epsilon + ) diff --git a/tests/eko/test_ev_operator.py b/tests/eko/test_ev_operator.py index 4a657e575..e242f5e17 100644 --- a/tests/eko/test_ev_operator.py +++ b/tests/eko/test_ev_operator.py @@ -6,7 +6,7 @@ from eko import anomalous_dimensions as ad from eko import basis_rotation as br from eko import interpolation, mellin -from eko.evolution_operator import Operator, gamma_ns_fact, gamma_singlet_fact, quad_ker +from eko.evolution_operator import Operator, quad_ker from eko.evolution_operator.grid import OperatorGrid from eko.interpolation import InterpolatorDispatcher from eko.kernels import non_singlet as ns @@ -15,32 +15,6 @@ from eko.thresholds import ThresholdsAtlas -def test_gamma_ns_fact(monkeypatch): - gamma_ns = np.array([1.0, 0.5, 0.25]) - monkeypatch.setattr(ad, "gamma_ns", lambda *args: gamma_ns.copy()) - gamma_ns_LO_0 = gamma_ns_fact(0, "NS_p", 1, 3, 0) - np.testing.assert_allclose(gamma_ns_LO_0, gamma_ns) - gamma_ns_LO_1 = gamma_ns_fact(0, "NS_p", 1, 3, 1) - np.testing.assert_allclose(gamma_ns_LO_1, gamma_ns) - gamma_ns_NLO_1 = gamma_ns_fact(1, "NS_p", 1, 3, 1) - assert gamma_ns_NLO_1[1] < gamma_ns[1] - gamma_ns_NNLO_1 = gamma_ns_fact(2, "NS_p", 1, 3, 1) - assert gamma_ns_NNLO_1[2] - gamma_ns[2] == 8.0 - - -def test_gamma_singlet_fact(monkeypatch): - gamma_s = np.array([1.0, 0.5, 0.25]) - monkeypatch.setattr(ad, "gamma_singlet", lambda *args: gamma_s.copy()) - gamma_s_LO_0 = gamma_singlet_fact(0, 1, 3, 0) - np.testing.assert_allclose(gamma_s_LO_0, gamma_s) - gamma_s_LO_1 = gamma_singlet_fact(0, 1, 3, 1) - np.testing.assert_allclose(gamma_s_LO_1, gamma_s) - gamma_s_NLO_1 = gamma_singlet_fact(1, 1, 3, 1) - assert gamma_s_NLO_1[1] < gamma_s[1] - gamma_s_NNLO_1 = gamma_singlet_fact(2, 1, 3, 1) - assert gamma_s_NNLO_1[2] - gamma_s[2] == 8.0 - - def test_quad_ker(monkeypatch): monkeypatch.setattr( mellin, "Talbot_path", lambda *args: 2 @@ -67,6 +41,7 @@ def test_quad_ker(monkeypatch): L=0, ev_op_iterations=0, ev_op_max_order=0, + sv_mode=0, ) np.testing.assert_allclose(res_ns, 0.0) res_s = quad_ker( @@ -83,6 +58,7 @@ def test_quad_ker(monkeypatch): L=0, ev_op_iterations=0, ev_op_max_order=0, + sv_mode=0, ) np.testing.assert_allclose(res_s, 1.0) res_s = quad_ker( @@ -99,8 +75,29 @@ def test_quad_ker(monkeypatch): L=0, ev_op_iterations=0, ev_op_max_order=0, + sv_mode=0, ) np.testing.assert_allclose(res_s, 0.0) + for mode in ["NS_p", "S_qq"]: + for sv in [1, 2]: + res_sv = quad_ker( + u=0, + order=0, + mode=mode, + method="", + is_log=True, + logx=1.0, + areas=np.zeros(3), + a1=1, + a0=2, + nf=3, + L=0, + ev_op_iterations=0, + ev_op_max_order=0, + sv_mode=sv, + ) + np.testing.assert_allclose(res_sv, 1.0) + monkeypatch.setattr(interpolation, "log_evaluate_Nx", lambda *args: 0) res_ns = quad_ker( u=0, @@ -116,6 +113,7 @@ def test_quad_ker(monkeypatch): L=0, ev_op_iterations=0, ev_op_max_order=0, + sv_mode=0, ) np.testing.assert_allclose(res_ns, 0.0) @@ -129,7 +127,7 @@ def test_labels(self): 1, 2, ) - assert sorted(o.labels()) == sorted( + assert sorted(o.labels) == sorted( ["NS_p", "NS_m", "NS_v", "S_qq", "S_qg", "S_gq", "S_gg"] ) o = Operator( @@ -139,7 +137,7 @@ def test_labels(self): 1, 2, ) - assert sorted(o.labels()) == [] + assert sorted(o.labels) == [] def test_compute(self, monkeypatch): # setup objs @@ -165,6 +163,7 @@ def test_compute(self, monkeypatch): "MaxNfPdf": 6, "MaxNfAs": 6, "HQ": "POLE", + "ModSV": None, } operators_card = { "Q2grid": [1, 10], @@ -218,13 +217,13 @@ def test_compute(self, monkeypatch): def test_pegasus_path(): def quad_ker_pegasus( - u, order, mode, method, logx, areas, a1, a0, nf, L, ev_op_iterations + u, order, mode, method, logx, areas, a1, a0, nf, ev_op_iterations ): # compute the mellin inversion as done in pegasus phi = 3 / 4 * np.pi c = 1.9 n = complex(c + u * np.exp(1j * phi)) - gamma_ns = gamma_ns_fact(order, mode, n, nf, L) + gamma_ns = ad.gamma_ns(order, mode[-1], n, nf) ker = ns.dispatcher( order, method, @@ -237,7 +236,7 @@ def quad_ker_pegasus( pj = interpolation.log_evaluate_Nx(n, logx, areas) return np.imag(np.exp(1j * phi) / np.pi * pj * ker) - # It might be useful to test with a different fuction + # It might be useful to test with a different function # monkeypatch.setattr(ns, "dispatcher", lambda x, *args: np.exp( - x ** 2 ) ) xgrid = np.geomspace(1e-7, 1, 10) int_disp = InterpolatorDispatcher(xgrid, 1, True) @@ -269,6 +268,7 @@ def quad_ker_pegasus( L, ev_op_iterations, 10, + 0, ), epsabs=1e-12, epsrel=1e-5, @@ -289,7 +289,6 @@ def quad_ker_pegasus( a1, a0, nf, - L, ev_op_iterations, ), epsabs=1e-12, diff --git a/tests/eko/test_msbar_masses.py b/tests/eko/test_msbar_masses.py index ec77fe302..beed1d40e 100644 --- a/tests/eko/test_msbar_masses.py +++ b/tests/eko/test_msbar_masses.py @@ -38,10 +38,10 @@ def test_compute_msbar_mass(self): for method in ["EXA", "EXP"]: for order in [1, 2, 3]: theory_dict.update({"ModEv": method, "PTO": order}) - strong_coupling = StrongCoupling.from_dict(theory_dict) # compute the scale such msbar(m) = m m2_computed = msbar_masses.compute(theory_dict) + strong_coupling = StrongCoupling.from_dict(theory_dict, m2_computed) m2_test = [] for nf in [3, 4, 5]: # compute msbar( m ) @@ -56,7 +56,7 @@ def test_compute_msbar_mass(self): q2_to=m2_computed[nf - 3], ) ) - np.testing.assert_allclose(m2_computed, m2_test, rtol=4e-3) + np.testing.assert_allclose(m2_computed, m2_test, rtol=6e-4) def test_compute_msbar_mass_VFNS(self): # test the solution now with some initial contition @@ -71,9 +71,9 @@ def test_compute_msbar_mass_VFNS(self): "Qmb": 85.0, } ) - strong_coupling = StrongCoupling.from_dict(theory_dict) # compute the scale such msbar(m) = m m2_computed = msbar_masses.compute(theory_dict) + strong_coupling = StrongCoupling.from_dict(theory_dict, m2_computed) m2_test = [] for nf in [3, 4, 5]: # compute msbar( m ) @@ -88,7 +88,7 @@ def test_compute_msbar_mass_VFNS(self): q2_to=m2_computed[nf - 3], ) ) - np.testing.assert_allclose(m2_computed, m2_test, rtol=3e-3) + np.testing.assert_allclose(m2_computed, m2_test, rtol=6e-4) def test_errors(self): diff --git a/tests/eko/test_ome.py b/tests/eko/test_ome.py index 9cd9c2105..dc8b6355c 100644 --- a/tests/eko/test_ome.py +++ b/tests/eko/test_ome.py @@ -150,7 +150,7 @@ def test_quad_ker(monkeypatch): else: np.testing.assert_allclose(res_ns, 0.0) - # test exact intrisic inverse kernel + # test exact intrinsic inverse kernel labels.extend( [ "S_Hq", @@ -214,12 +214,13 @@ class TestOperatorMatrixElement: "mc": 1.0, "mb": 4.75, "mt": 173.0, - "kcThr": np.inf, - "kbThr": np.inf, + "kcThr": 1.0, + "kbThr": 1.0, "ktThr": np.inf, "MaxNfPdf": 6, "MaxNfAs": 6, "HQ": "POLE", + "ModSV": None, } def test_labels(self): @@ -278,7 +279,7 @@ def test_compute_lo(self): InterpolatorDispatcher.from_dict(operators_card), ) o = OperatorMatrixElement(g.config, g.managers, is_backward=False) - o.compute(self.theory_card["mb"] ** 2, L=0, is_msbar=False) + o.compute(self.theory_card["mb"] ** 2, nf=4, L=0, is_msbar=False) dim = o.ome_members["NS_qq"].value.shape for idx in ["S", "NS"]: @@ -317,12 +318,12 @@ def test_compute_nlo(self): g = OperatorGrid.from_dict( t, operators_card, - ThresholdsAtlas.from_dict(self.theory_card), - StrongCoupling.from_dict(self.theory_card), + ThresholdsAtlas.from_dict(t), + StrongCoupling.from_dict(t), InterpolatorDispatcher.from_dict(operators_card), ) o = OperatorMatrixElement(g.config, g.managers, is_backward=False) - o.compute(self.theory_card["mb"] ** 2, L=0, is_msbar=False) + o.compute(t["mb"] ** 2, nf=4, L=0, is_msbar=False) dim = len(operators_card["interpolation_xgrid"]) shape = (dim, dim) diff --git a/tests/eko/test_runner.py b/tests/eko/test_runner.py index b31f88331..1523e4034 100644 --- a/tests/eko/test_runner.py +++ b/tests/eko/test_runner.py @@ -30,6 +30,7 @@ "Qmc": 1.0, "Qmb": 4.75, "Qmt": 173.0, + "ModSV": None, } operators_card = { "Q2grid": [10, 100], diff --git a/tests/eko/test_strong_coupling.py b/tests/eko/test_strong_coupling.py index 4b6e1b88c..26751138d 100644 --- a/tests/eko/test_strong_coupling.py +++ b/tests/eko/test_strong_coupling.py @@ -28,6 +28,7 @@ def test_from_dict(self): "ktThr": 1.0, "MaxNfAs": 6, "HQ": "POLE", + "ModSV": None, } sc = StrongCoupling.from_dict(d) assert sc.a_s(d["Qref"] ** 2) == d["alphas"] / (4.0 * np.pi) diff --git a/tests/eko/test_sv_expanded.py b/tests/eko/test_sv_expanded.py new file mode 100644 index 000000000..076df98f0 --- /dev/null +++ b/tests/eko/test_sv_expanded.py @@ -0,0 +1,99 @@ +# -*- coding: utf-8 -*- + +import numpy as np + +from eko.anomalous_dimensions import gamma_ns, gamma_singlet +from eko.beta import beta_0 +from eko.kernels import non_singlet, singlet +from eko.scale_variations import expanded, exponentiated + + +def test_ns_sv_dispacher(): + """Test to identity""" + order = 3 + gamma_ns = np.random.rand(order + 1) + L = 0 + nf = 5 + a_s = 0.35 + np.testing.assert_allclose( + expanded.non_singlet_variation(gamma_ns, a_s, order, nf, L), 1.0 + ) + + +def test_singlet_sv_dispacher(): + """Test to identity""" + order = 3 + gamma_singlet = np.random.rand(order + 1, 2, 2) + L = 0 + nf = 5 + a_s = 0.35 + np.testing.assert_allclose( + expanded.singlet_variation(gamma_singlet, a_s, order, nf, L), np.eye(2) + ) + + +def test_scale_variation_a_vs_b(): + r""" + Test ``ModSV=exponentiated`` kernel vs ``ModSV=expanded``. + We test that the quantity :math:`(ker_A - ker_B)/ker_{unv}` depends + only on the accuracy in the :math:`\alpha_s` expansion + and not in the size of the `fact_to_ren` itself. + However in our implementation the exponentiated mode depends on + the actual value of a0 and a1, since the evolution integral in :math:`\alpha_s` + is evaluated. Thus this test ratio :math:`(ker_A - ker_B)/ker_{unv}` + still contains a dependency on `fact_to_ren`. + """ + nf = 5 + n = 10 + a1 = 0.118 / (4 * np.pi) + a0 = 0.2 / (4 * np.pi) + method = "truncated" + + def scheme_diff(g, k, pto, is_singlet): + """ + :math:`(ker_A - ker_B)/ker_{unv}` for truncated expansion + Effects due to non commutativity are neglected thus, + he accuracy of singlet quantities is slightly worst. + """ + if pto >= 1: + diff = g[0] * k * a0 + if pto >= 2: + b0 = beta_0(nf) + g02 = g[0] @ g[0] if is_singlet else g[0] ** 2 + diff += a0**2 * g[1] * k - k**2 * ( + 1 / 2 * a0**2 * b0 * g[0] + a1 * a0 * g02 - 1 / 2 * a0**2 * g02 + ) + return diff + + # TODO: perform this test also at N3LO, once evolution kernels + # will be implemented + + for L in [np.log(0.5), np.log(2)]: + for order in [1, 2]: + # Non singlet kernels + gns = gamma_ns(order, "p", n, nf) + ker = non_singlet.dispatcher( + order, method, gns, a1, a0, nf, ev_op_iterations=1 + ) + gns_a = exponentiated.gamma_variation(gns.copy(), order, nf, L) + ker_a = non_singlet.dispatcher( + order, method, gns_a, a1, a0, nf, ev_op_iterations=1 + ) + ker_b = ker * expanded.non_singlet_variation(gns, a1, order, nf, L) + ns_diff = scheme_diff(gns, L, order, False) + np.testing.assert_allclose((ker_a - ker_b) / ker, ns_diff, atol=1e-3) + + # Singlet kernels + gs = gamma_singlet(order, n, nf) + ker = singlet.dispatcher( + order, method, gs, a1, a0, nf, ev_op_iterations=1, ev_op_max_order=1 + ) + gs_a = exponentiated.gamma_variation(gs.copy(), order, nf, L) + ker_a = singlet.dispatcher( + order, method, gs_a, a1, a0, nf, ev_op_iterations=1, ev_op_max_order=1 + ) + ker_b = ker @ expanded.singlet_variation(gs, a1, order, nf, L) + s_diff = scheme_diff(gs, L, order, True) + np.testing.assert_allclose( + (ker_a - ker_b) @ np.linalg.inv(ker), s_diff, atol=5e-3 + ) diff --git a/tests/eko/test_sv_exponentiated.py b/tests/eko/test_sv_exponentiated.py new file mode 100644 index 000000000..fe0367825 --- /dev/null +++ b/tests/eko/test_sv_exponentiated.py @@ -0,0 +1,33 @@ +# -*- coding: utf-8 -*- + +import numpy as np + +from eko.scale_variations.exponentiated import gamma_variation + + +def test_gamma_ns_fact(): + gamma_ns = np.array([1.0, 0.5, 0.25, 0.125]) + gamma_ns_LO_0 = gamma_variation(gamma_ns.copy(), 0, 3, 0) + np.testing.assert_allclose(gamma_ns_LO_0, gamma_ns) + gamma_ns_LO_1 = gamma_variation(gamma_ns.copy(), 0, 3, 1) + np.testing.assert_allclose(gamma_ns_LO_1, gamma_ns) + gamma_ns_NLO_1 = gamma_variation(gamma_ns.copy(), 1, 3, 1) + assert gamma_ns_NLO_1[1] < gamma_ns[1] + gamma_ns_NNLO_1 = gamma_variation(gamma_ns.copy(), 2, 3, 1) + assert gamma_ns_NNLO_1[2] - gamma_ns[2] == 8.0 + gamma_ns_N3LO_0 = gamma_variation(gamma_ns.copy(), 3, 3, 0) + assert gamma_ns_N3LO_0[3] == gamma_ns[3] + + +def test_gamma_singlet_fact(): + gamma_s = np.array([1.0, 0.5, 0.25, 0.125]) + gamma_s_LO_0 = gamma_variation(gamma_s.copy(), 0, 3, 0) + np.testing.assert_allclose(gamma_s_LO_0, gamma_s) + gamma_s_LO_1 = gamma_variation(gamma_s.copy(), 0, 3, 1) + np.testing.assert_allclose(gamma_s_LO_1, gamma_s) + gamma_s_NLO_1 = gamma_variation(gamma_s.copy(), 1, 3, 1) + assert gamma_s_NLO_1[1] < gamma_s[1] + gamma_s_NNLO_1 = gamma_variation(gamma_s.copy(), 2, 3, 1) + assert gamma_s_NNLO_1[2] - gamma_s[2] == 8.0 + gamma_s_N3LO_0 = gamma_variation(gamma_s.copy(), 3, 3, 0) + assert gamma_s_N3LO_0[3] == gamma_s[3] diff --git a/tests/eko/test_thresholds.py b/tests/eko/test_thresholds.py index b5b55ff7a..99eefb312 100644 --- a/tests/eko/test_thresholds.py +++ b/tests/eko/test_thresholds.py @@ -5,7 +5,7 @@ import numpy as np import pytest -from eko.thresholds import PathSegment, ThresholdsAtlas +from eko.thresholds import PathSegment, ThresholdsAtlas, flavor_shift, is_downward_path class TestPathSegment: @@ -197,3 +197,23 @@ def test_errors(self): HQ="FAIL", ), ) + + +def test_downward_path(): + thr_atlas = ThresholdsAtlas( + masses=np.power([2, 3, 4], 2), + q2_ref=91**2, + nf_ref=3, + thresholds_ratios=[1, 1, 1], + ) + q2_to = 5**2 + path_3 = thr_atlas.path(q2_to, nf_to=3) + # path_3 is downward in q2 + is_downward = is_downward_path(path_3) + assert is_downward is True + assert flavor_shift(is_downward) == 4 + # path_6 is downward in q2, but forward in nf + path_6 = thr_atlas.path(q2_to, nf_to=6) + is_downward = is_downward_path(path_6) + assert is_downward is False + assert flavor_shift(is_downward) == 3