From 78e533cee78e30e441973fdc219d4be8f27ac4eb Mon Sep 17 00:00:00 2001 From: juacrumar Date: Tue, 16 Jul 2024 09:41:34 +0200 Subject: [PATCH 1/3] use nf to generate alpha_s --- src/ekobox/info_file.py | 24 +++++++++--------------- 1 file changed, 9 insertions(+), 15 deletions(-) diff --git a/src/ekobox/info_file.py b/src/ekobox/info_file.py index 26e3cb170..fb744ce4b 100644 --- a/src/ekobox/info_file.py +++ b/src/ekobox/info_file.py @@ -68,19 +68,13 @@ def build( hqm_scheme=theory_card.heavy.masses_scheme, thresholds_ratios=np.power(list(iter(theory_card.heavy.matching_ratios)), 2.0), ) - alphas_values = np.array( - [ - 4.0 - * np.pi - * sc.a_s( - muf2, - ) - for muf2 in operators_card.mu2grid - ], - dtype=float, - ) - template_info["AlphaS_Vals"] = alphas_values.tolist() - template_info["AlphaS_Qs"] = np.array( - [mu for mu, _ in operators_card.mugrid] - ).tolist() + + alphas_values = [] + alphas_qs = [] + for mu, nf in operators_card.mugrid: + alphas_values.append(float(4.0 * np.pi * sc.a_s(mu * mu, nf_to=nf))) + alphas_qs.append(mu) + + template_info["AlphaS_Vals"] = alphas_values + template_info["AlphaS_Qs"] = alphas_qs return template_info From cb6953df9adeef319c6f23f6aedcc21993ad601b Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 16 Jul 2024 14:15:58 +0300 Subject: [PATCH 2/3] Check sorting of alphas in LHAPDF info --- src/ekobox/evol_pdf.py | 10 +----- src/ekobox/info_file.py | 62 +++++++++++++++++++++++++++------- src/ekobox/utils.py | 9 +++++ tests/ekobox/test_info_file.py | 40 ++++++++++++++++++++++ 4 files changed, 100 insertions(+), 21 deletions(-) diff --git a/src/ekobox/evol_pdf.py b/src/ekobox/evol_pdf.py index ddd80f75a..0045cdd28 100644 --- a/src/ekobox/evol_pdf.py +++ b/src/ekobox/evol_pdf.py @@ -1,13 +1,13 @@ """Tools to evolve actual PDFs.""" import pathlib -from collections import defaultdict from eko import basis_rotation as br from eko.io import EKO from eko.runner import managed from . import apply, genpdf, info_file +from .utils import regroup_evolgrid DEFAULT_NAME = "eko.tar" @@ -95,14 +95,6 @@ def evolve_pdfs( genpdf.install_pdf(name) -def regroup_evolgrid(evolgrid: list): - """Split evolution points by nf and sort by scale.""" - by_nf = defaultdict(list) - for q, nf in sorted(evolgrid, key=lambda ep: ep[1]): - by_nf[nf].append(q) - return {nf: sorted(qs) for nf, qs in by_nf.items()} - - def collect_blocks(evolved_PDF: dict, q2block_per_nf: dict, xgrid: list): """Collect all LHAPDF blocks for a given replica. diff --git a/src/ekobox/info_file.py b/src/ekobox/info_file.py index fb744ce4b..c9ac88b0c 100644 --- a/src/ekobox/info_file.py +++ b/src/ekobox/info_file.py @@ -9,6 +9,7 @@ from eko.io.runcards import OperatorCard, TheoryCard from .genpdf import load +from .utils import regroup_evolgrid def build( @@ -16,19 +17,19 @@ def build( operators_card: OperatorCard, num_members: int, info_update: dict, -): - """Generate a lhapdf info file from theory and operators card. +) -> dict: + """Generate a lhapdf info file. Parameters ---------- - theory_card : dict + theory_card : theory card - operators_card : dict - operators_card - num_members : int + operators_card : + operators card + num_members : number of pdf set members - info_update : dict - info to update + info_update : + additional info to update Returns ------- @@ -56,8 +57,45 @@ def build( template_info["MCharm"] = theory_card.heavy.masses.c.value template_info["MBottom"] = theory_card.heavy.masses.b.value template_info["MTop"] = theory_card.heavy.masses.t.value + # dump alphas + template_info.update(build_alphas(theory_card, operators_card)) + return template_info + + +def build_alphas( + theory_card: TheoryCard, + operators_card: OperatorCard, +) -> dict: + """Generate a couplings section of lhapdf info file. + + Parameters + ---------- + theory_card : dict + theory card + operators_card : dict + operators card + + Returns + ------- + dict + info file section in lhapdf format + """ + # start with meta stuff + template_info = {} template_info["AlphaS_MZ"] = theory_card.couplings.alphas template_info["AlphaS_OrderQCD"] = theory_card.order[0] - 1 + + # check we have disjoint scale ranges + evolgrid = regroup_evolgrid(operators_card.mugrid) + nfs = list(evolgrid.keys()) + for j in range(len(nfs) - 1): + # equal points are allowed by LHAPDF + if evolgrid[nfs[j]][-1] > evolgrid[nfs[j + 1]][0]: + raise ValueError( + f"Last scale point for nf={nfs[j]} is bigger than first in nf={nfs[j+1]}" + ) + + # add actual values evmod = couplings.couplings_mod_ev(operators_card.configs.evolution_method) quark_masses = [(x.value) ** 2 for x in theory_card.heavy.masses] sc = couplings.Couplings( @@ -68,12 +106,12 @@ def build( hqm_scheme=theory_card.heavy.masses_scheme, thresholds_ratios=np.power(list(iter(theory_card.heavy.matching_ratios)), 2.0), ) - alphas_values = [] alphas_qs = [] - for mu, nf in operators_card.mugrid: - alphas_values.append(float(4.0 * np.pi * sc.a_s(mu * mu, nf_to=nf))) - alphas_qs.append(mu) + for nf, mus in evolgrid.items(): + for mu in mus: + alphas_values.append(float(4.0 * np.pi * sc.a_s(mu * mu, nf_to=nf))) + alphas_qs.append(mu) template_info["AlphaS_Vals"] = alphas_values template_info["AlphaS_Qs"] = alphas_qs diff --git a/src/ekobox/utils.py b/src/ekobox/utils.py index 0047d0872..091cde15e 100644 --- a/src/ekobox/utils.py +++ b/src/ekobox/utils.py @@ -1,6 +1,7 @@ """Generic utilities to work with EKOs.""" import os +from collections import defaultdict from typing import Optional import numpy as np @@ -78,3 +79,11 @@ def ekos_product( if path is not None: final_eko.close() + + +def regroup_evolgrid(evolgrid: list): + """Split evolution points by nf and sort by scale.""" + by_nf = defaultdict(list) + for q, nf in sorted(evolgrid, key=lambda ep: ep[1]): + by_nf[nf].append(q) + return {nf: sorted(qs) for nf, qs in by_nf.items()} diff --git a/tests/ekobox/test_info_file.py b/tests/ekobox/test_info_file.py index 64af8eb93..f7267fbd4 100644 --- a/tests/ekobox/test_info_file.py +++ b/tests/ekobox/test_info_file.py @@ -1,6 +1,7 @@ import math import numpy as np +import pytest from ekobox import cards, info_file @@ -23,3 +24,42 @@ def test_build(): np.testing.assert_allclose(info["QMin"], math.sqrt(op.mu2grid[0]), rtol=1e-5) assert info["XMin"] == op.xgrid.raw[0] assert info["XMax"] == op.xgrid.raw[-1] == 1.0 + + +def test_build_alphas_good(): + """Good configurations.""" + theory = cards.example.theory() + theory.order = (2, 0) + theory.couplings.alphas = 0.2 + op = cards.example.operator() + op.mu0 = 1.0 + # base case + op.mugrid = [(10.0, 5), (100.0, 5)] + info = info_file.build_alphas(theory, op) + assert len(info["AlphaS_Vals"]) == 2 + np.testing.assert_allclose(info["AlphaS_Qs"], [10.0, 100.0]) + # several nf + op.mugrid = [(1.0, 3), (5.0, 3), (5.0, 4), (10.0, 5), (10.0, 5), (100.0, 5)] + info = info_file.build_alphas(theory, op) + assert len(info["AlphaS_Vals"]) == 6 + np.testing.assert_allclose(info["AlphaS_Qs"], [1.0, 5.0, 5.0, 10.0, 10.0, 100.0]) + # several nf with gap + op.mugrid = [(1.0, 3), (10.0, 3), (10.0, 5), (100.0, 5)] + info = info_file.build_alphas(theory, op) + assert len(info["AlphaS_Vals"]) == 4 + np.testing.assert_allclose(info["AlphaS_Qs"], [1.0, 10.0, 10.0, 100.0]) + + +def test_build_alphas_bad(): + """Bad configurations.""" + theory = cards.example.theory() + theory.order = (2, 0) + theory.couplings.alphas = 0.2 + op = cards.example.operator() + op.mu0 = 1.0 + op.mugrid = [(5.0, 3), (15.0, 4), (10.0, 5), (100.0, 5)] + with pytest.raises(ValueError, match="is bigger"): + info_file.build_alphas(theory, op) + op.mugrid = [(5.0, 3), (10.0, 3), (10.0, 4), (15.0, 4), (10.0, 5)] + with pytest.raises(ValueError, match="is bigger"): + info_file.build_alphas(theory, op) From d5f63fe7ac995fd9d0d38dba48f2a7969de01c28 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 16 Jul 2024 14:53:40 +0300 Subject: [PATCH 3/3] Lift scale check for LHAPDF --- src/ekobox/evol_pdf.py | 19 ++++++++++++++++--- src/ekobox/info_file.py | 13 ++----------- tests/ekobox/test_evol_pdf.py | 16 ++++++++++++++++ tests/ekobox/test_info_file.py | 20 ++------------------ 4 files changed, 36 insertions(+), 32 deletions(-) diff --git a/src/ekobox/evol_pdf.py b/src/ekobox/evol_pdf.py index 0045cdd28..abdb2ec6a 100644 --- a/src/ekobox/evol_pdf.py +++ b/src/ekobox/evol_pdf.py @@ -2,6 +2,8 @@ import pathlib +import numpy as np + from eko import basis_rotation as br from eko.io import EKO from eko.runner import managed @@ -46,6 +48,18 @@ def evolve_pdfs( info_update : dict dict of info to add or update to default info file """ + # separate by nf the evolgrid (and order per nf/q) + q2block_per_nf = regroup_evolgrid(operators_card.mugrid) + + # check we have disjoint scale ranges + nfs = list(q2block_per_nf.keys()) + for j in range(len(nfs) - 1): + # equal points are allowed by LHAPDF + if q2block_per_nf[nfs[j]][-1] > q2block_per_nf[nfs[j + 1]][0]: + raise ValueError( + f"Last scale point for nf={nfs[j]} is bigger than first in nf={nfs[j+1]}" + ) + # update op and th cards if path is not None: eko_path = pathlib.Path(path) @@ -79,9 +93,8 @@ def evolve_pdfs( info_update=info_update, ) - # separate by nf the evolgrid (and order per nf/q) - q2block_per_nf = regroup_evolgrid(eko_output.evolgrid) - + # in the eko scales are squared + q2block_per_nf = {nf: np.power(q2s, 2) for nf, q2s in q2block_per_nf.items()} # write all replicas all_member_blocks = [] for evolved_PDF in evolved_PDF_list: diff --git a/src/ekobox/info_file.py b/src/ekobox/info_file.py index c9ac88b0c..fc5c17f1d 100644 --- a/src/ekobox/info_file.py +++ b/src/ekobox/info_file.py @@ -84,18 +84,8 @@ def build_alphas( template_info = {} template_info["AlphaS_MZ"] = theory_card.couplings.alphas template_info["AlphaS_OrderQCD"] = theory_card.order[0] - 1 - - # check we have disjoint scale ranges + # prepare evolgrid = regroup_evolgrid(operators_card.mugrid) - nfs = list(evolgrid.keys()) - for j in range(len(nfs) - 1): - # equal points are allowed by LHAPDF - if evolgrid[nfs[j]][-1] > evolgrid[nfs[j + 1]][0]: - raise ValueError( - f"Last scale point for nf={nfs[j]} is bigger than first in nf={nfs[j+1]}" - ) - - # add actual values evmod = couplings.couplings_mod_ev(operators_card.configs.evolution_method) quark_masses = [(x.value) ** 2 for x in theory_card.heavy.masses] sc = couplings.Couplings( @@ -106,6 +96,7 @@ def build_alphas( hqm_scheme=theory_card.heavy.masses_scheme, thresholds_ratios=np.power(list(iter(theory_card.heavy.matching_ratios)), 2.0), ) + # add actual values alphas_values = [] alphas_qs = [] for nf, mus in evolgrid.items(): diff --git a/tests/ekobox/test_evol_pdf.py b/tests/ekobox/test_evol_pdf.py index 54bce6623..141787b49 100644 --- a/tests/ekobox/test_evol_pdf.py +++ b/tests/ekobox/test_evol_pdf.py @@ -1,4 +1,5 @@ import numpy as np +import pytest from banana import toy import eko @@ -50,6 +51,21 @@ def test_evolve_pdfs_dump_path(fake_lhapdf, cd): assert p.exists() +def test_evolve_pdfs_bad_scales(fake_lhapdf, cd): + """Bad scale configurations.""" + theory, op = init_cards() + n = "test_evolve_pdfs_bad_scales" + op = cards.example.operator() + op.mugrid = [(5.0, 3), (15.0, 4), (10.0, 5), (100.0, 5)] + with pytest.raises(ValueError, match="is bigger"): + with cd(fake_lhapdf): + ev_p.evolve_pdfs([toy.mkPDF("", 0)], theory, op, name=n, path=fake_lhapdf) + op.mugrid = [(5.0, 3), (10.0, 3), (10.0, 4), (15.0, 4), (10.0, 5)] + with pytest.raises(ValueError, match="is bigger"): + with cd(fake_lhapdf): + ev_p.evolve_pdfs([toy.mkPDF("", 0)], theory, op, name=n, path=fake_lhapdf) + + def test_evolve_pdfs_dump_file(fake_lhapdf, cd): theory, op = init_cards() n = "test_evolve_pdfs_dump_file" diff --git a/tests/ekobox/test_info_file.py b/tests/ekobox/test_info_file.py index f7267fbd4..acf3f7272 100644 --- a/tests/ekobox/test_info_file.py +++ b/tests/ekobox/test_info_file.py @@ -1,7 +1,6 @@ import math import numpy as np -import pytest from ekobox import cards, info_file @@ -34,12 +33,12 @@ def test_build_alphas_good(): op = cards.example.operator() op.mu0 = 1.0 # base case - op.mugrid = [(10.0, 5), (100.0, 5)] + op.mugrid = [(100.0, 5), (10.0, 5)] info = info_file.build_alphas(theory, op) assert len(info["AlphaS_Vals"]) == 2 np.testing.assert_allclose(info["AlphaS_Qs"], [10.0, 100.0]) # several nf - op.mugrid = [(1.0, 3), (5.0, 3), (5.0, 4), (10.0, 5), (10.0, 5), (100.0, 5)] + op.mugrid = [(5.0, 4), (10.0, 5), (1.0, 3), (5.0, 3), (10.0, 5), (100.0, 5)] info = info_file.build_alphas(theory, op) assert len(info["AlphaS_Vals"]) == 6 np.testing.assert_allclose(info["AlphaS_Qs"], [1.0, 5.0, 5.0, 10.0, 10.0, 100.0]) @@ -48,18 +47,3 @@ def test_build_alphas_good(): info = info_file.build_alphas(theory, op) assert len(info["AlphaS_Vals"]) == 4 np.testing.assert_allclose(info["AlphaS_Qs"], [1.0, 10.0, 10.0, 100.0]) - - -def test_build_alphas_bad(): - """Bad configurations.""" - theory = cards.example.theory() - theory.order = (2, 0) - theory.couplings.alphas = 0.2 - op = cards.example.operator() - op.mu0 = 1.0 - op.mugrid = [(5.0, 3), (15.0, 4), (10.0, 5), (100.0, 5)] - with pytest.raises(ValueError, match="is bigger"): - info_file.build_alphas(theory, op) - op.mugrid = [(5.0, 3), (10.0, 3), (10.0, 4), (15.0, 4), (10.0, 5)] - with pytest.raises(ValueError, match="is bigger"): - info_file.build_alphas(theory, op)