Skip to content

Commit

Permalink
allow convolution with 2 ekos
Browse files Browse the repository at this point in the history
  • Loading branch information
@giacomomagni
giacomomagni committed Jun 3, 2024
1 parent e91b2cd commit 57dff37
Show file tree
Hide file tree
Showing 2 changed files with 110 additions and 41 deletions.
59 changes: 39 additions & 20 deletions src/pineko/evolve.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -180,13 +180,13 @@ def write_operator_card(pineappl_grid, default_card, card_path, tcard):
)

# For hardonic obs we might need to dump 2 eko cards
operators_card["configs"]["polarized"] = eko1 == "polPDF"
if eko1 == eko2:
with open(card_path, "w", encoding="UTF-8") as f:
yaml.safe_dump(operators_card, f)
f.write(f"# {pineko_version=}")
else:
# dump eko1
operators_card["configs"]["polarized"] = eko1 == "polPDF"
card_name = f"{card_path.stem}_eko1.yaml"
with open(card_path.parent / card_name, "w", encoding="UTF-8") as f:
yaml.safe_dump(operators_card, f)
Expand All @@ -203,12 +203,13 @@ def write_operator_card(pineappl_grid, default_card, card_path, tcard):

def evolve_grid(
grid,
operators,
operators_a,
fktable_path,
max_as,
max_al,
xir,
xif,
operators_b=None,
assumptions="Nf6Ind",
comparison_pdf=None,
meta_data=None,
Expand All @@ -220,7 +221,7 @@ def evolve_grid(
----------
grid : pineappl.grid.Grid
unconvoluted grid
operators : eko.EKO
operators_a : eko.EKO
evolution operator
fktable_path : str
target path for convoluted grid
Expand All @@ -232,6 +233,8 @@ def evolve_grid(
renormalization scale variation
xif : float
factorization scale variation
operators_b: eko.EKO
additonal evolution operator if different from operators_a
assumptions : str
assumptions on the flavor dimension
comparison_pdf : None or str
Expand All @@ -253,22 +256,29 @@ def evolve_grid(
evol_info = grid.evolve_info(order_mask)
x_grid = evol_info.x1
mur2_grid = evol_info.ren1
xif = 1.0 if operators.operator_card.configs.scvar_method is not None else xif
tcard = operators.theory_card
opcard = operators.operator_card
xif = 1.0 if operators_a.operator_card.configs.scvar_method is not None else xif
tcard = operators_a.theory_card
opcard = operators_a.operator_card
# rotate the targetgrid
if "integrability_version" in grid.key_values():
x_grid = np.append(x_grid, 1.0)
eko.io.manipulate.xgrid_reshape(
operators, targetgrid=eko.interpolation.XGrid(x_grid)
)
check.check_grid_and_eko_compatible(grid, operators, xif, max_as, max_al)
# rotate to evolution (if doable and necessary)
if np.allclose(operators.bases.inputpids, br.flavor_basis_pids):
eko.io.manipulate.to_evol(operators)
# Here we are checking if the EKO contains the rotation matrix (flavor to evol)
elif not np.allclose(operators.bases.inputpids, br.rotate_flavor_to_evolution):
raise ValueError("The EKO is neither in flavor nor in evolution basis.")

def xgrid_reshape(operators):
eko.io.manipulate.xgrid_reshape(
operators, targetgrid=eko.interpolation.XGrid(x_grid)
)
check.check_grid_and_eko_compatible(grid, operators, xif, max_as, max_al)
# rotate to evolution (if doable and necessary)
if np.allclose(operators.bases.inputpids, br.flavor_basis_pids):
eko.io.manipulate.to_evol(operators)
# Here we are checking if the EKO contains the rotation matrix (flavor to evol)
elif not np.allclose(operators.bases.inputpids, br.rotate_flavor_to_evolution):
raise ValueError("The EKO is neither in flavor nor in evolution basis.")

xgrid_reshape(operators_a)
if operators_b is not None:
xgrid_reshape(operators_b)

# PineAPPL wants alpha_s = 4*pi*a_s
# remember that we already accounted for xif in the opcard generation
evmod = eko.couplings.couplings_mod_ev(opcard.configs.evolution_method)
Expand Down Expand Up @@ -296,18 +306,27 @@ def evolve_grid(
]
# We need to use ekompatibility in order to pass a dictionary to pineappl
fktable = grid.evolve(
ekompatibility.pineappl_layout(operators),
ekompatibility.pineappl_layout(operators_a),
xir * xir * mur2_grid,
alphas_values,
"evol",
operators_b=ekompatibility.pineappl_layout(operators_b),
order_mask=order_mask,
xi=(xir, xif),
)
rich.print(f"Optimizing for {assumptions}")
fktable.optimize(assumptions)
fktable.set_key_value("eko_version", operators.metadata.version)
fktable.set_key_value("eko_theory_card", json.dumps(operators.theory_card.raw))
fktable.set_key_value("eko_operator_card", json.dumps(operators.operator_card.raw))
fktable.set_key_value("eko_version", operators_a.metadata.version)
fktable.set_key_value("eko_theory_card", json.dumps(operators_a.theory_card.raw))

fktable.set_key_value(
"eko_operator_card", json.dumps(operators_a.operator_card.raw)
)
if operators_b is not None:
fktable.set_key_value(
"eko_operator_card_b", json.dumps(operators_b.operator_card.raw)
)

fktable.set_key_value("pineko_version", version.__version__)
if meta_data is not None:
for k, v in meta_data.items():
Expand Down
92 changes: 71 additions & 21 deletions src/pineko/theory.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -398,8 +398,27 @@ def fk(self, name, grid_path, tcard, pdf):
grid = pineappl.grid.Grid.read(grid_path)
# remove zero subgrid
grid.optimize()

# Do you need one or multiple ekos?
kv = grid.key_values()
if "convolution_type_1" in kv:
eko1 = kv["convolution_type_1"]
# TODO: this case is now deprecated and should be remved from yadism and pinefarm
elif "polarized" in kv:
eko1 = "polPDF"
else:
eko1 = "PDF"
eko2 = kv.get("convolution_type_2", "PDF")

# setup data
eko_filename = self.ekos_path() / f"{name}.tar"
if eko1 == eko2:
eko_filename = [self.ekos_path() / f"{name}.tar"]
else:
eko_filename = [
self.ekos_path() / f"{name}_eko1.tar",
self.ekos_path() / f"{name}_eko2.tar",
]

fk_filename = self.fks_path / f"{name}.{parser.EXT}"
if fk_filename.exists():
if not self.overwrite:
Expand Down Expand Up @@ -434,16 +453,28 @@ def fk(self, name, grid_path, tcard, pdf):
if not np.isclose(xif, 1.0):
check_scvar_evolve(grid, max_as, max_al, check.Scale.FACT)
# loading ekos to produce a tmp copy
with eko.EKO.read(eko_filename) as operators:
n_ekos = len(eko_filename)
with eko.EKO.read(eko_filename[0]) as operators_a:

# Skip the computation of the fktable if the eko is empty
if len(operators.mu2grid) == 0 and check.is_num_fonll(tcard["FNS"]):
if len(operators_a.mu2grid) == 0 and check.is_num_fonll(tcard["FNS"]):
rich.print("[green] Skipping empty eko for nFONLL.")
return
eko_tmp_path = (
operators.paths.root.parent / f"eko-tmp-{name}-{np.random.rand()}.tar"

eko_tmp_path_a = (
operators_a.paths.root.parent / f"eko-tmp-{name}-{np.random.rand()}.tar"
)
operators.deepcopy(eko_tmp_path)
with eko.EKO.edit(eko_tmp_path) as operators:
operators_a.deepcopy(eko_tmp_path_a)

if n_ekos > 1:
with eko.EKO.read(eko_filename[1]) as operators_b:
eko_tmp_path_b = (
operators_a.paths.root.parent
/ f"eko-tmp-{name}-{np.random.rand()}.tar"
)
operators_b.deepcopy(eko_tmp_path_b)

with eko.EKO.edit(eko_tmp_path_a) as operators_a:
# Obtain the assumptions hash
assumptions = theory_card.construct_assumptions(tcard)
# do it!
Expand All @@ -456,7 +487,6 @@ def fk(self, name, grid_path, tcard, pdf):
xif,
)
start_time = time.perf_counter()

rich.print(
rich.panel.Panel.fit(
"Computing ...", style="magenta", box=rich.box.SQUARE
Expand All @@ -466,20 +496,40 @@ def fk(self, name, grid_path, tcard, pdf):
f"= {fk_filename}\n",
f"with max_as={max_as}, max_al={max_al}, xir={xir}, xif={xif}",
)
_grid, _fk, comparison = evolve.evolve_grid(
grid,
operators,
fk_filename,
max_as,
max_al,
xir=xir,
xif=xif,
assumptions=assumptions,
comparison_pdf=pdf,
meta_data={"theory_card": json.dumps(tcard)},
)

if n_ekos == 1:
_grid, _fk, comparison = evolve.evolve_grid(
grid,
operators_a,
fk_filename,
max_as,
max_al,
xir=xir,
xif=xif,
assumptions=assumptions,
comparison_pdf=pdf,
meta_data={"theory_card": json.dumps(tcard)},
)
else:
with eko.EKO.edit(eko_tmp_path_b) as operators_b:
_grid, _fk, comparison = evolve.evolve_grid(
grid,
operators_a,
fk_filename,
max_as,
max_al,
xir=xir,
xif=xif,
assumptions=assumptions,
operators_b=operators_b,
comparison_pdf=pdf,
meta_data={"theory_card": json.dumps(tcard)},
)
# Remove tmp ekos
eko_tmp_path_b.unlink()

# Remove tmp ekos
eko_tmp_path.unlink()
eko_tmp_path_a.unlink()

logger.info(
"Finished computation of %s - took %f s",
Expand Down

0 comments on commit 57dff37

Please sign in to comment.