From 557a5fb4a209304f233d6879878739edab04e0b4 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 11 Jul 2022 11:32:10 +0200 Subject: [PATCH 01/24] Branch debug stuff --- debug.py | 1 + 1 file changed, 1 insertion(+) diff --git a/debug.py b/debug.py index 9e506f05..73c92414 100644 --- a/debug.py +++ b/debug.py @@ -12,6 +12,7 @@ bo = eko.output.Output.load_tar("data/ekos/2201/test.tar") bpdf = apply.apply_pdf(bo, gonly) +# a is the output of the sv script plot(np.log(a[:, 2]), a[:, 4]) plot(np.log(a[:, 2]), a[:, 6]) plot(np.log(a[:, 2]), 4 * 5 / 9 * a[:, 2] * bpdf[300.0]["pdfs"][2]) From 9a116bdf8969fb5005fc1bafa482f9ea281a4528 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Thu, 8 Sep 2022 18:15:51 +0200 Subject: [PATCH 02/24] Reconcile pre-commit --- debug2.py | 65 ++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 43 insertions(+), 22 deletions(-) diff --git a/debug2.py b/debug2.py index 9b4617b4..2f76c0f0 100644 --- a/debug2.py +++ b/debug2.py @@ -1,52 +1,73 @@ # -*- coding: utf-8 -*- +import numpy as np import yadism import yaml -import numpy as np -from scipy.integrate import quad -from scipy.interpolate import interp1d from eko.anomalous_dimensions import as1 from eko.harmonics import S1 from eko.interpolation import InterpolatorDispatcher +# from scipy.integrate import quad +# from scipy.interpolate import interp1d # def mel(y, n): # xs = out["interpolation_xgrid"] # f = interp1d(xs, y) # return quad(lambda x: x ** (n - 1.0) * f(x), min(xs), 1.0) + def mel2(out, ys, n): interp = InterpolatorDispatcher.from_dict(out) lnxmin = np.log(out["interpolation_xgrid"][0]) - res = 0. + res = 0.0 for y, bf in zip(ys, interp): pj = bf(n, lnxmin) * np.exp(n * lnxmin) res += y * pj return res -with open("210.yaml","r") as f: t = yaml.safe_load(f) -with open("test.yaml","r") as f: o = yaml.safe_load(f) -#out = yadism.output.Output.load_tar("test.tar") -out = yadism.run_yadism(t,o) +with open("210.yaml", "r") as f: + t = yaml.safe_load(f) +with open("test.yaml", "r") as f: + o = yaml.safe_load(f) +# out = yadism.output.Output.load_tar("test.tar") +out = yadism.run_yadism(t, o) out.dump_tar("test.tar") -locbf10 = [out["F2_light"][j].orders[(0,0,0,0)][0][-3][10] for j in range(len(out["interpolation_xgrid"]))] -pqqcbf10 = [out["F2_light"][j].orders[(1,0,0,1)][0][-3][10] for j in range(len(out["interpolation_xgrid"]))] -pgqgbf10 = [out["F2_light"][j].orders[(1,0,0,1)][0][7][10] for j in range(len(out["interpolation_xgrid"]))] +locbf10 = [ + out["F2_light"][j].orders[(0, 0, 0, 0)][0][-3][10] + for j in range(len(out["interpolation_xgrid"])) +] +pqqcbf10 = [ + out["F2_light"][j].orders[(1, 0, 0, 1)][0][-3][10] + for j in range(len(out["interpolation_xgrid"])) +] +pgqgbf10 = [ + out["F2_light"][j].orders[(1, 0, 0, 1)][0][7][10] + for j in range(len(out["interpolation_xgrid"])) +] -mel2bf10 = np.array([mel2(out,[0]*10 + [1] + [0]*39,n) for n in range(1,10)]) +mel2bf10 = np.array([mel2(out, [0] * 10 + [1] + [0] * 39, n) for n in range(1, 10)]) -mel2lobf10 = np.array([mel2(out,locbf10,n) for n in range(1,10)]) -mel2pqqcbf10 = np.array([mel2(out,pqqcbf10,n) for n in range(1,10)]) -mel2pgqgbf10 = np.array([mel2(out,pgqgbf10,n) for n in range(1,10)]) +mel2lobf10 = np.array([mel2(out, locbf10, n) for n in range(1, 10)]) +mel2pqqcbf10 = np.array([mel2(out, pqqcbf10, n) for n in range(1, 10)]) +mel2pgqgbf10 = np.array([mel2(out, pgqgbf10, n) for n in range(1, 10)]) -locbf20 = [out["F2_light"][j].orders[(0,0,0,0)][0][-3][20] for j in range(len(out["interpolation_xgrid"]))] -pqqcbf20 = [out["F2_light"][j].orders[(1,0,0,1)][0][-3][20] for j in range(len(out["interpolation_xgrid"]))] -pgqgbf20 = [out["F2_light"][j].orders[(1,0,0,1)][0][7][20] for j in range(len(out["interpolation_xgrid"]))] +locbf20 = [ + out["F2_light"][j].orders[(0, 0, 0, 0)][0][-3][20] + for j in range(len(out["interpolation_xgrid"])) +] +pqqcbf20 = [ + out["F2_light"][j].orders[(1, 0, 0, 1)][0][-3][20] + for j in range(len(out["interpolation_xgrid"])) +] +pgqgbf20 = [ + out["F2_light"][j].orders[(1, 0, 0, 1)][0][7][20] + for j in range(len(out["interpolation_xgrid"])) +] -mel2bf20 = np.array([mel2(out,[0]*20 + [1] + [0]*29,n) for n in range(1,10)]) -mel2lobf20 = np.array([mel2(out,locbf20,n) for n in range(1,10)]) -mel2pqqcbf20 = np.array([mel2(out,pqqcbf20,n) for n in range(1,10)]) -mel2pgqgbf20 = np.array([mel2(out,pgqgbf20,n) for n in range(1,10)]) +mel2bf20 = np.array([mel2(out, [0] * 20 + [1] + [0] * 29, n) for n in range(1, 10)]) +mel2lobf20 = np.array([mel2(out, locbf20, n) for n in range(1, 10)]) +mel2pqqcbf20 = np.array([mel2(out, pqqcbf20, n) for n in range(1, 10)]) +mel2pgqgbf20 = np.array([mel2(out, pgqgbf20, n) for n in range(1, 10)]) gns = np.array([as1.gamma_ns(n, S1(n)) for n in range(1, 10)]) ggq = np.array([np.nan if n == 1 else as1.gamma_gq(n) for n in range(1, 10)]) From 91258c6b5fb5838fec6315f71230774d93a466b5 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 21 Sep 2022 18:43:25 +0200 Subject: [PATCH 03/24] Add yet more debug stuff --- debug3.py | 83 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 debug3.py diff --git a/debug3.py b/debug3.py new file mode 100644 index 00000000..1259220d --- /dev/null +++ b/debug3.py @@ -0,0 +1,83 @@ +# -*- coding: utf-8 -*- +import lhapdf +import numpy as np +import pineappl +import yadism +from eko import interpolation +from yadism.coefficient_functions.splitting_functions import lo, nlo +from yadism.esf.conv import convolute_operator + +mode_g = not True + +if mode_g: + pdf = lhapdf.mkPDF("gonly", 0) +else: + pdf = lhapdf.mkPDF("conly", 0) + + +b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/test.pineappl.lz4") +c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/test.pineappl.lz4") +out = yadism.output.Output.load_tar("test0.tar") + +bb1 = b1.convolute_with_one(2212, pdf.xfxQ2) +cc1 = c1.convolute_with_one(2212, pdf.xfxQ2) + +loc = np.array( + [ + out["F2_light"][j].orders[(0, 0, 0, 0)][0][-3] + for j in range(len(out["interpolation_xgrid"])) + ] +) +nloc = np.array( + [ + out["F2_light"][j].orders[(1, 0, 0, 0)][0][-3] + for j in range(len(out["interpolation_xgrid"])) + ] +) +nlog = np.array( + [ + out["F2_light"][j].orders[(1, 0, 0, 0)][0][7] + for j in range(len(out["interpolation_xgrid"])) + ] +) +pqqc = np.array( + [ + out["F2_light"][j].orders[(1, 0, 0, 1)][0][-3] + for j in range(len(out["interpolation_xgrid"])) + ] +) +pqgg = np.array( + [ + out["F2_light"][j].orders[(1, 0, 0, 1)][0][7] + for j in range(len(out["interpolation_xgrid"])) + ] +) +interp = interpolation.InterpolatorDispatcher.from_dict(out, mode_N=False) +pqq = convolute_operator(lo.pqq(5), interp)[0] +pqg = convolute_operator(lo.pqg(5), interp)[0] +pgq = convolute_operator(nlo.pgq0(5), interp)[0] +pgg = convolute_operator(nlo.pgg0(5), interp)[0] + +muf2 = 300.0 * 2.0**2 +# a = pdf.alphasQ2(muf2) +a = 0.0078 +if mode_g: + f = np.array([pdf.xfxQ2(21, x, muf2) / x for x in out["interpolation_xgrid"]]) +else: + f = np.array([pdf.xfxQ2(4, x, muf2) / x for x in out["interpolation_xgrid"]]) + +# if mode_g: +# # this yields 4/9 as it should +# print((pqgg @ f)/(pqg.T @ f)/out["interpolation_xgrid"]) +# else: +# # this yields 22/90 ?! +# print((pqqc @ f)/(pqq.T @ f)/out["interpolation_xgrid"]) +# print(cc1/((loc + a*(nloc + pqqc)) @ f)) + +diff_grid = bb1 - cc1 +if mode_g: + diff_ana = a**2 * np.log(2.0**2) * ((nloc @ pqg.T @ f) + (nlog @ pgg.T @ f)) +else: + diff_ana = a**2 * np.log(2.0**2) * ((nloc @ pqq.T @ f) + (nlog @ pgq.T @ f)) + +print(diff_grid / diff_ana) From 61339cff775bebf94dcb2684b02e9684a7e1bc31 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Thu, 22 Sep 2022 13:01:21 +0200 Subject: [PATCH 04/24] Add more debug --- debug3.py | 44 ++++++++++++++++++++++++++++++++++---------- 1 file changed, 34 insertions(+), 10 deletions(-) diff --git a/debug3.py b/debug3.py index 1259220d..9377b314 100644 --- a/debug3.py +++ b/debug3.py @@ -1,22 +1,27 @@ # -*- coding: utf-8 -*- import lhapdf import numpy as np +import matplotlib.pyplot as plt import pineappl import yadism from eko import interpolation from yadism.coefficient_functions.splitting_functions import lo, nlo from yadism.esf.conv import convolute_operator -mode_g = not True +mode_g = True +mode_no_nlo = not True if mode_g: pdf = lhapdf.mkPDF("gonly", 0) else: pdf = lhapdf.mkPDF("conly", 0) - -b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/test.pineappl.lz4") -c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/test.pineappl.lz4") +if mode_no_nlo: + b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/no-nlo-test.pineappl.lz4") + c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/no-nlo-test.pineappl.lz4") +else: + b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/test.pineappl.lz4") + c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/test.pineappl.lz4") out = yadism.output.Output.load_tar("test0.tar") bb1 = b1.convolute_with_one(2212, pdf.xfxQ2) @@ -60,7 +65,12 @@ muf2 = 300.0 * 2.0**2 # a = pdf.alphasQ2(muf2) -a = 0.0078 +a = 0.011 + +if mode_no_nlo: + nlog *= 0 + nloc *= 0 + if mode_g: f = np.array([pdf.xfxQ2(21, x, muf2) / x for x in out["interpolation_xgrid"]]) else: @@ -70,14 +80,28 @@ # # this yields 4/9 as it should # print((pqgg @ f)/(pqg.T @ f)/out["interpolation_xgrid"]) # else: -# # this yields 22/90 ?! +# # this yields 22/90 = 2*(1+4+1+4+1)/9/(2*5) as it should # print((pqqc @ f)/(pqq.T @ f)/out["interpolation_xgrid"]) -# print(cc1/((loc + a*(nloc + pqqc)) @ f)) + +# if mode_g: +# print(cc1/((a*(nlog + pqgg*np.log(1./2.**2))) @ f)) +# else: +# print(cc1/((loc + a*(nloc + pqqc*np.log(1./2.**2))) @ f)) diff_grid = bb1 - cc1 +# plt.plot(diff_grid, label="grid") +# nloc *= 1.9 +# nlog *= 1.14 if mode_g: - diff_ana = a**2 * np.log(2.0**2) * ((nloc @ pqg.T @ f) + (nlog @ pgg.T @ f)) + # plt.plot(a**2 * np.log(1./2.0**2) * ((nloc @ pqg.T @ f)), label="pqg") + # plt.plot(a**2 * np.log(1./2.0**2) * ((nlog @ pgg.T @ f)), label="pgg") + diff_ana = a**2 * np.log(1./2.0**2) * ((nloc @ pqg.T @ f) + (nlog @ pgg.T @ f)) else: - diff_ana = a**2 * np.log(2.0**2) * ((nloc @ pqq.T @ f) + (nlog @ pgq.T @ f)) - + # plt.plot(a**2 * np.log(1./2.0**2) * ((nloc @ pqq.T @ f)), label="pqq") + # plt.plot(a**2 * np.log(1./2.0**2) * ((nlog @ pgq.T @ f)), label="pgq") + diff_ana = a**2 * np.log(1./2.0**2) * ((nloc @ pqq.T @ f) + (nlog @ pgq.T @ f)) +# plt.legend() +# plt.show() +# print(diff_grid/cc1) +# print(diff_grid) print(diff_grid / diff_ana) From 0def85acaa1c0794c98608a865562974a465db68 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Thu, 22 Sep 2022 18:19:30 +0200 Subject: [PATCH 05/24] Happy debugging --- debug3.py | 86 +++++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 74 insertions(+), 12 deletions(-) diff --git a/debug3.py b/debug3.py index 9377b314..352b3804 100644 --- a/debug3.py +++ b/debug3.py @@ -1,21 +1,36 @@ # -*- coding: utf-8 -*- +import sys import lhapdf import numpy as np import matplotlib.pyplot as plt import pineappl import yadism from eko import interpolation +from eko import output +from ekomark.apply import apply_pdf + +# from eko.output.struct import EKO from yadism.coefficient_functions.splitting_functions import lo, nlo from yadism.esf.conv import convolute_operator +# mode can be gluon only or charm only mode_g = True -mode_no_nlo = not True +if len(sys.argv) > 1: + fl = sys.argv[1].strip() + if fl == "g": + mode_g = True + if fl == "c": + mode_g = False + +mode_no_nlo = True +# load respective PDF set if mode_g: pdf = lhapdf.mkPDF("gonly", 0) else: pdf = lhapdf.mkPDF("conly", 0) +# load FK tables if mode_no_nlo: b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/no-nlo-test.pineappl.lz4") c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/no-nlo-test.pineappl.lz4") @@ -24,9 +39,11 @@ c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/test.pineappl.lz4") out = yadism.output.Output.load_tar("test0.tar") +# compute predictions from FK bb1 = b1.convolute_with_one(2212, pdf.xfxQ2) cc1 = c1.convolute_with_one(2212, pdf.xfxQ2) +# load ingredients from yadism loc = np.array( [ out["F2_light"][j].orders[(0, 0, 0, 0)][0][-3] @@ -57,51 +74,96 @@ for j in range(len(out["interpolation_xgrid"])) ] ) + +# compute splitting functions interp = interpolation.InterpolatorDispatcher.from_dict(out, mode_N=False) pqq = convolute_operator(lo.pqq(5), interp)[0] -pqg = convolute_operator(lo.pqg(5), interp)[0] +pqg = convolute_operator(lo.pqg(5), interp)[0] # / 10 # * 55 / 100 pgq = convolute_operator(nlo.pgq0(5), interp)[0] pgg = convolute_operator(nlo.pgg0(5), interp)[0] muf2 = 300.0 * 2.0**2 # a = pdf.alphasQ2(muf2) a = 0.011 +L = np.log(1.0 / 2.0**2) if mode_no_nlo: nlog *= 0 nloc *= 0 +# compute PDFs if mode_g: f = np.array([pdf.xfxQ2(21, x, muf2) / x for x in out["interpolation_xgrid"]]) else: f = np.array([pdf.xfxQ2(4, x, muf2) / x for x in out["interpolation_xgrid"]]) +# check yadism SV is proportional splitting functions +# print("check yadism SV") # if mode_g: # # this yields 4/9 as it should # print((pqgg @ f)/(pqg.T @ f)/out["interpolation_xgrid"]) # else: -# # this yields 22/90 = 2*(1+4+1+4+1)/9/(2*5) as it should +# # this yields 22/90 = 2*(1+4+1+4+1)/9/(2*5) # print((pqqc @ f)/(pqq.T @ f)/out["interpolation_xgrid"]) +# extract a by comparing C against yadism +# print("check C result") # if mode_g: -# print(cc1/((a*(nlog + pqgg*np.log(1./2.**2))) @ f)) +# print(cc1/((a*(nlog + pqgg*L)) @ f)) # else: -# print(cc1/((loc + a*(nloc + pqqc*np.log(1./2.**2))) @ f)) +# print(cc1/((loc + a*(nloc + pqqc*L)) @ f)) + +# construct K +Kqq = np.eye(len(out["interpolation_xgrid"])) + a * L * pqq.T +Kgg = np.eye(len(out["interpolation_xgrid"])) + a * L * pgg.T +Kqg = a * L * pqg.T +Kgq = a * L * pgq.T + +# check B EKO is K +print("check B EKO result") +ekob = output.Output.load_tar("data/ekos/2205/test.tar") +fb = apply_pdf(ekob, pdf) +if mode_g: + gb = Kgg @ f + cb = Kqg @ f +else: + gb = Kgq @ f + cb = Kqq @ f +print(fb[300.0]["pdfs"][21] / gb) +print(fb[300.0]["pdfs"][4] / cb) + +# check C EKO is identity +# print("check C EKO result") +# ekoc = output.Output.load_tar("data/ekos/3205/test.tar") +# fc = apply_pdf(ekoc, pdf) +# if mode_g: +# print(fc[2.**2 * 300.]["pdfs"][21]/f) +# print(fc[2.**2 * 300.]["pdfs"][4] - np.zeros_like(f)) +# else: +# print(fc[2.**2 * 300.]["pdfs"][4]/f) +# print(fc[2.**2 * 300.]["pdfs"][21] - np.zeros_like(f)) + +# check B result +print("check B result") +if mode_g: + print(bb1 / (((loc + a * nloc) @ Kqg + (a * nlog) @ Kgg) @ f)) +else: + print(bb1 / (((loc + a * nloc) @ Kqq + (a * nlog) @ Kgq) @ f)) diff_grid = bb1 - cc1 # plt.plot(diff_grid, label="grid") # nloc *= 1.9 # nlog *= 1.14 if mode_g: - # plt.plot(a**2 * np.log(1./2.0**2) * ((nloc @ pqg.T @ f)), label="pqg") - # plt.plot(a**2 * np.log(1./2.0**2) * ((nlog @ pgg.T @ f)), label="pgg") - diff_ana = a**2 * np.log(1./2.0**2) * ((nloc @ pqg.T @ f) + (nlog @ pgg.T @ f)) + # plt.plot(a**2 * L * ((nloc @ pqg.T @ f)), label="pqg") + # plt.plot(a**2 * L * ((nlog @ pgg.T @ f)), label="pgg") + diff_ana = a**2 * L * ((nloc @ pqg.T @ f) + (nlog @ pgg.T @ f)) else: - # plt.plot(a**2 * np.log(1./2.0**2) * ((nloc @ pqq.T @ f)), label="pqq") - # plt.plot(a**2 * np.log(1./2.0**2) * ((nlog @ pgq.T @ f)), label="pgq") - diff_ana = a**2 * np.log(1./2.0**2) * ((nloc @ pqq.T @ f) + (nlog @ pgq.T @ f)) + # plt.plot(a**2 * L * ((nloc @ pqq.T @ f)), label="pqq") + # plt.plot(a**2 * L * ((nlog @ pgq.T @ f)), label="pgq") + diff_ana = a**2 * L * ((nloc @ pqq.T @ f) + (nlog @ pgq.T @ f)) # plt.legend() # plt.show() # print(diff_grid/cc1) # print(diff_grid) -print(diff_grid / diff_ana) +# print(diff_grid / diff_ana) From deaf045ffaf2eab5b5e9c1759adeced9a23644c6 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Fri, 23 Sep 2022 18:45:06 +0200 Subject: [PATCH 06/24] Consider all quarks in debug --- debug3.py | 98 +++++++++++++++++++++++++++++++++---------------------- 1 file changed, 59 insertions(+), 39 deletions(-) diff --git a/debug3.py b/debug3.py index 352b3804..725397e7 100644 --- a/debug3.py +++ b/debug3.py @@ -23,6 +23,8 @@ mode_g = False mode_no_nlo = True +if len(sys.argv) > 2: + mode_no_nlo = sys.argv[1].strip().lower() == "nlo" # load respective PDF set if mode_g: @@ -50,12 +52,24 @@ for j in range(len(out["interpolation_xgrid"])) ] ) +los = np.array( + [ + out["F2_light"][j].orders[(0, 0, 0, 0)][0][-4] + for j in range(len(out["interpolation_xgrid"])) + ] +) nloc = np.array( [ out["F2_light"][j].orders[(1, 0, 0, 0)][0][-3] for j in range(len(out["interpolation_xgrid"])) ] ) +nlos = np.array( + [ + out["F2_light"][j].orders[(1, 0, 0, 0)][0][-3] + for j in range(len(out["interpolation_xgrid"])) + ] +) nlog = np.array( [ out["F2_light"][j].orders[(1, 0, 0, 0)][0][7] @@ -78,7 +92,7 @@ # compute splitting functions interp = interpolation.InterpolatorDispatcher.from_dict(out, mode_N=False) pqq = convolute_operator(lo.pqq(5), interp)[0] -pqg = convolute_operator(lo.pqg(5), interp)[0] # / 10 # * 55 / 100 +pqg = convolute_operator(lo.pqg(5), interp)[0] / 10 # * 55 / 100 pgq = convolute_operator(nlo.pgq0(5), interp)[0] pgg = convolute_operator(nlo.pgg0(5), interp)[0] @@ -87,9 +101,21 @@ a = 0.011 L = np.log(1.0 / 2.0**2) +# shut yadism references down if mode_no_nlo: nlog *= 0 nloc *= 0 + nlos *= 0 + +# create joint quark terms +loq = 2 * (3 * los + 2 * loc) +nloq = 2 * (3 * nlos + 2 * nloc) + +# construct K +Kqq = np.eye(len(out["interpolation_xgrid"])) + a * L * pqq.T +Kgg = np.eye(len(out["interpolation_xgrid"])) + a * L * pgg.T +Kqg = a * L * pqg.T +Kgq = a * L * pgq.T # compute PDFs if mode_g: @@ -100,11 +126,9 @@ # check yadism SV is proportional splitting functions # print("check yadism SV") # if mode_g: -# # this yields 4/9 as it should -# print((pqgg @ f)/(pqg.T @ f)/out["interpolation_xgrid"]) +# print((pqgg @ f)/(loq @ pqg.T @ f)) # else: -# # this yields 22/90 = 2*(1+4+1+4+1)/9/(2*5) -# print((pqqc @ f)/(pqq.T @ f)/out["interpolation_xgrid"]) +# print((pqqc @ f)/(loc @ pqq.T @ f)) # extract a by comparing C against yadism # print("check C result") @@ -113,24 +137,18 @@ # else: # print(cc1/((loc + a*(nloc + pqqc*L)) @ f)) -# construct K -Kqq = np.eye(len(out["interpolation_xgrid"])) + a * L * pqq.T -Kgg = np.eye(len(out["interpolation_xgrid"])) + a * L * pgg.T -Kqg = a * L * pqg.T -Kgq = a * L * pgq.T - # check B EKO is K -print("check B EKO result") -ekob = output.Output.load_tar("data/ekos/2205/test.tar") -fb = apply_pdf(ekob, pdf) -if mode_g: - gb = Kgg @ f - cb = Kqg @ f -else: - gb = Kgq @ f - cb = Kqq @ f -print(fb[300.0]["pdfs"][21] / gb) -print(fb[300.0]["pdfs"][4] / cb) +# print("check B EKO result") +# ekob = output.Output.load_tar("data/ekos/2205/test.tar") +# fb = apply_pdf(ekob, pdf) +# if mode_g: +# gb = Kgg @ f +# cb = Kqg @ f +# else: +# gb = Kgq @ f +# cb = Kqq @ f +# print(fb[300.0]["pdfs"][21] / gb) +# print(fb[300.0]["pdfs"][4] / cb) # check C EKO is identity # print("check C EKO result") @@ -146,24 +164,26 @@ # check B result print("check B result") if mode_g: - print(bb1 / (((loc + a * nloc) @ Kqg + (a * nlog) @ Kgg) @ f)) + print(bb1 / (((loq + a * nloq) @ Kqg + (a * nlog) @ Kgg) @ f)) else: print(bb1 / (((loc + a * nloc) @ Kqq + (a * nlog) @ Kgq) @ f)) -diff_grid = bb1 - cc1 -# plt.plot(diff_grid, label="grid") -# nloc *= 1.9 -# nlog *= 1.14 -if mode_g: - # plt.plot(a**2 * L * ((nloc @ pqg.T @ f)), label="pqg") - # plt.plot(a**2 * L * ((nlog @ pgg.T @ f)), label="pgg") - diff_ana = a**2 * L * ((nloc @ pqg.T @ f) + (nlog @ pgg.T @ f)) -else: - # plt.plot(a**2 * L * ((nloc @ pqq.T @ f)), label="pqq") - # plt.plot(a**2 * L * ((nlog @ pgq.T @ f)), label="pgq") - diff_ana = a**2 * L * ((nloc @ pqq.T @ f) + (nlog @ pgq.T @ f)) -# plt.legend() -# plt.show() -# print(diff_grid/cc1) -# print(diff_grid) +# check difference between B and C +# print("check B - C") +# diff_grid = bb1 - cc1 +# # plt.plot(diff_grid, label="grid") +# # nloc *= 1.9 +# # nlog *= 1.14 +# if mode_g: +# # plt.plot(a**2 * L * ((nloc @ pqg.T @ f)), label="pqg") +# # plt.plot(a**2 * L * ((nlog @ pgg.T @ f)), label="pgg") +# diff_ana = a**2 * L * ((nloq @ pqg.T @ f) + (nlog @ pgg.T @ f)) +# else: +# # plt.plot(a**2 * L * ((nloc @ pqq.T @ f)), label="pqq") +# # plt.plot(a**2 * L * ((nlog @ pgq.T @ f)), label="pgq") +# diff_ana = a**2 * L * ((nloc @ pqq.T @ f) + (nlog @ pgq.T @ f)) +# # plt.legend() +# # plt.show() +# # print(diff_grid/cc1) +# # print(diff_grid) # print(diff_grid / diff_ana) From 7e7abacf9744c7117bc3d8325da43eccf8039a5f Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 5 Oct 2022 14:40:42 +0200 Subject: [PATCH 07/24] Predict B-C --- debug3.py | 81 +++++++++++++++++++++++++++++++++---------------------- 1 file changed, 49 insertions(+), 32 deletions(-) diff --git a/debug3.py b/debug3.py index 725397e7..6fe7b66a 100644 --- a/debug3.py +++ b/debug3.py @@ -5,8 +5,11 @@ import matplotlib.pyplot as plt import pineappl import yadism +import yaml +import eko from eko import interpolation from eko import output +from eko import strong_coupling as sc from ekomark.apply import apply_pdf # from eko.output.struct import EKO @@ -96,9 +99,15 @@ pgq = convolute_operator(nlo.pgq0(5), interp)[0] pgg = convolute_operator(nlo.pgg0(5), interp)[0] -muf2 = 300.0 * 2.0**2 -# a = pdf.alphasQ2(muf2) -a = 0.011 +q2 = 300.0 +muf2 = q2 * 2.0**2 +# load theory card +with open("data/theory_cards/2205.yaml") as fd: + tc = yaml.safe_load(fd) +astrong = sc.StrongCoupling.from_dict(tc) +# actually, there are two alpha_s values in the game! +ac = astrong.a_s(muf2) +ab = astrong.a_s(q2) L = np.log(1.0 / 2.0**2) # shut yadism references down @@ -112,10 +121,10 @@ nloq = 2 * (3 * nlos + 2 * nloc) # construct K -Kqq = np.eye(len(out["interpolation_xgrid"])) + a * L * pqq.T -Kgg = np.eye(len(out["interpolation_xgrid"])) + a * L * pgg.T -Kqg = a * L * pqg.T -Kgq = a * L * pgq.T +Kqq = np.eye(len(out["interpolation_xgrid"])) + ac * L * pqq.T +Kgg = np.eye(len(out["interpolation_xgrid"])) + ac * L * pgg.T +Kqg = ac * L * pqg.T +Kgq = ac * L * pgq.T # compute PDFs if mode_g: @@ -133,9 +142,9 @@ # extract a by comparing C against yadism # print("check C result") # if mode_g: -# print(cc1/((a*(nlog + pqgg*L)) @ f)) +# print(cc1/((ac*(nlog + pqgg*L)) @ f)) # else: -# print(cc1/((loc + a*(nloc + pqqc*L)) @ f)) +# print(cc1/((loc + ac*(nloc + pqqc*L)) @ f)) # check B EKO is K # print("check B EKO result") @@ -162,28 +171,36 @@ # print(fc[2.**2 * 300.]["pdfs"][21] - np.zeros_like(f)) # check B result -print("check B result") -if mode_g: - print(bb1 / (((loq + a * nloq) @ Kqg + (a * nlog) @ Kgg) @ f)) -else: - print(bb1 / (((loc + a * nloc) @ Kqq + (a * nlog) @ Kgq) @ f)) - -# check difference between B and C -# print("check B - C") -# diff_grid = bb1 - cc1 -# # plt.plot(diff_grid, label="grid") -# # nloc *= 1.9 -# # nlog *= 1.14 +# print("check B result") # if mode_g: -# # plt.plot(a**2 * L * ((nloc @ pqg.T @ f)), label="pqg") -# # plt.plot(a**2 * L * ((nlog @ pgg.T @ f)), label="pgg") -# diff_ana = a**2 * L * ((nloq @ pqg.T @ f) + (nlog @ pgg.T @ f)) +# print(bb1 / (((loq + ab * nloq) @ Kqg + (ab * nlog) @ Kgg) @ f)) +# plt.plot(bb1) +# plt.plot((((loq + ab * nloq) @ Kqg + (ab * nlog) @ Kgg) @ f)) # else: -# # plt.plot(a**2 * L * ((nloc @ pqq.T @ f)), label="pqq") -# # plt.plot(a**2 * L * ((nlog @ pgq.T @ f)), label="pgq") -# diff_ana = a**2 * L * ((nloc @ pqq.T @ f) + (nlog @ pgq.T @ f)) -# # plt.legend() -# # plt.show() -# # print(diff_grid/cc1) -# # print(diff_grid) -# print(diff_grid / diff_ana) +# print(bb1 / (((loc + ab * nloc) @ Kqq + (ab * nlog) @ Kgq) @ f)) +# plt.show() + +# check difference between B and C +print("check B - C") +diff_grid = bb1 - cc1 +# plt.plot(bb1, label="b-grid") +# plt.plot(cc1, label="c-grid") +# nloc *= 1.9 +# nlog *= 1.14 +if mode_g: + # plt.plot(a**2 * L * ((nloc @ pqg.T @ f)), label="pqg") + # plt.plot(a**2 * L * ((nlog @ pgg.T @ f)), label="pgg") + diff_ana = (nlog * (ab - ac) + ac * ab * L * ((nloq @ pqg.T) + (nlog @ pgg.T))) @ f +else: + # plt.plot(a**2 * L * ((nloc @ pqq.T @ f)), label="pqq") + # plt.plot(a**2 * L * ((nlog @ pgq.T @ f)), label="pgq") + diff_ana = (nloc * (ab - ac) + ab * ac * L * ((nloc @ pqq.T) + (nlog @ pgq.T))) @ f +# print(diff_grid/cc1) +print(bb1 / cc1) +print(diff_grid / diff_ana) +plt.plot(diff_grid, label="grid") +plt.plot(diff_ana, label="ana") +# plt.plot((diff_grid/bb1-1)*100, label="rel err. diff to b") +# plt.plot(diff_grid/diff_ana, label="ratio") +plt.legend() +plt.show() From df5c4b60fcf7529f8e64589bc8b02040234294a5 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Thu, 6 Oct 2022 17:53:37 +0200 Subject: [PATCH 08/24] Start NNLO debug --- debug3.py | 217 +++++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 166 insertions(+), 51 deletions(-) diff --git a/debug3.py b/debug3.py index 6fe7b66a..f08377f0 100644 --- a/debug3.py +++ b/debug3.py @@ -1,4 +1,5 @@ # -*- coding: utf-8 -*- +from enum import Enum import sys import lhapdf import numpy as np @@ -11,6 +12,7 @@ from eko import output from eko import strong_coupling as sc from ekomark.apply import apply_pdf +from eko.beta import beta_0 # from eko.output.struct import EKO from yadism.coefficient_functions.splitting_functions import lo, nlo @@ -19,15 +21,26 @@ # mode can be gluon only or charm only mode_g = True if len(sys.argv) > 1: - fl = sys.argv[1].strip() + fl = sys.argv[1].strip().lower() if fl == "g": mode_g = True if fl == "c": mode_g = False -mode_no_nlo = True +# order mode +class Order(Enum): + nonlo = 0 + nlo = 1 + nnlo = 2 + + +order_mode = Order.nonlo if len(sys.argv) > 2: - mode_no_nlo = sys.argv[1].strip().lower() == "nlo" + om = sys.argv[1].strip().lower() + if om == "nlo": + order_mode = Order.nlo + elif om == "nnlo": + order_mode = Order.nnlo # load respective PDF set if mode_g: @@ -36,19 +49,23 @@ pdf = lhapdf.mkPDF("conly", 0) # load FK tables -if mode_no_nlo: +if order_mode == Order.nonlo: b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/no-nlo-test.pineappl.lz4") c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/no-nlo-test.pineappl.lz4") -else: +elif order_mode == Order.nlo: b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/test.pineappl.lz4") c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/test.pineappl.lz4") +elif order_mode == Order.nnlo: + b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/test-nnlo.pineappl.lz4") + c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/test-nnlo.pineappl.lz4") out = yadism.output.Output.load_tar("test0.tar") # compute predictions from FK bb1 = b1.convolute_with_one(2212, pdf.xfxQ2) cc1 = c1.convolute_with_one(2212, pdf.xfxQ2) -# load ingredients from yadism +# load ingredients from yadism --- +# LO loc = np.array( [ out["F2_light"][j].orders[(0, 0, 0, 0)][0][-3] @@ -61,6 +78,7 @@ for j in range(len(out["interpolation_xgrid"])) ] ) +# NLO nloc = np.array( [ out["F2_light"][j].orders[(1, 0, 0, 0)][0][-3] @@ -79,25 +97,73 @@ for j in range(len(out["interpolation_xgrid"])) ] ) -pqqc = np.array( +# NLO SV +nlosvc = np.array( [ out["F2_light"][j].orders[(1, 0, 0, 1)][0][-3] for j in range(len(out["interpolation_xgrid"])) ] ) -pqgg = np.array( +nlosvg = np.array( [ out["F2_light"][j].orders[(1, 0, 0, 1)][0][7] for j in range(len(out["interpolation_xgrid"])) ] ) +# NNLO +nnloc = np.array( + [ + out["F2_light"][j].orders[(2, 0, 0, 0)][0][-3] + for j in range(len(out["interpolation_xgrid"])) + ] +) +nnlos = np.array( + [ + out["F2_light"][j].orders[(2, 0, 0, 0)][0][-3] + for j in range(len(out["interpolation_xgrid"])) + ] +) +nnlog = np.array( + [ + out["F2_light"][j].orders[(2, 0, 0, 0)][0][7] + for j in range(len(out["interpolation_xgrid"])) + ] +) +# NNLO SV +nnlosv1c = np.array( + [ + out["F2_light"][j].orders[(2, 0, 0, 1)][0][-3] + for j in range(len(out["interpolation_xgrid"])) + ] +) +nnlosv1g = np.array( + [ + out["F2_light"][j].orders[(2, 0, 0, 1)][0][7] + for j in range(len(out["interpolation_xgrid"])) + ] +) +nnlosv2c = np.array( + [ + out["F2_light"][j].orders[(2, 0, 0, 2)][0][-3] + for j in range(len(out["interpolation_xgrid"])) + ] +) +nnlosv2g = np.array( + [ + out["F2_light"][j].orders[(2, 0, 0, 2)][0][7] + for j in range(len(out["interpolation_xgrid"])) + ] +) # compute splitting functions interp = interpolation.InterpolatorDispatcher.from_dict(out, mode_N=False) -pqq = convolute_operator(lo.pqq(5), interp)[0] -pqg = convolute_operator(lo.pqg(5), interp)[0] / 10 # * 55 / 100 -pgq = convolute_operator(nlo.pgq0(5), interp)[0] -pgg = convolute_operator(nlo.pgg0(5), interp)[0] +pqq0 = convolute_operator(lo.pqq(5), interp)[0].T +pqg0 = convolute_operator(lo.pqg(5), interp)[0].T / 10 # * 55 / 100 +pgq0 = convolute_operator(nlo.pgq0(5), interp)[0].T +pgg0 = convolute_operator(nlo.pgg0(5), interp)[0].T +pqq1 = convolute_operator(nlo.pqq1(5), interp)[0].T +pqg1 = convolute_operator(nlo.pqg1(5), interp)[0].T +beta0I = beta_0(5) * np.eye(pqg0.shape[0]) q2 = 300.0 muf2 = q2 * 2.0**2 @@ -111,20 +177,25 @@ L = np.log(1.0 / 2.0**2) # shut yadism references down -if mode_no_nlo: +if order_mode == Order.nonlo: nlog *= 0 nloc *= 0 nlos *= 0 +if order_mode == Order.nonlo or order_mode == Order.nlo: + nnlog *= 0 + nnloc *= 0 + nnlos *= 0 # create joint quark terms loq = 2 * (3 * los + 2 * loc) nloq = 2 * (3 * nlos + 2 * nloc) +nnloq = 2 * (3 * nnlos + 2 * nnloc) # construct K -Kqq = np.eye(len(out["interpolation_xgrid"])) + ac * L * pqq.T -Kgg = np.eye(len(out["interpolation_xgrid"])) + ac * L * pgg.T -Kqg = ac * L * pqg.T -Kgq = ac * L * pgq.T +Kqq = np.eye(len(out["interpolation_xgrid"])) + ac * L * pqq0 +Kgg = np.eye(len(out["interpolation_xgrid"])) + ac * L * pgg0 +Kqg = ac * L * pqg0 +Kgq = ac * L * pgq0 # compute PDFs if mode_g: @@ -132,19 +203,64 @@ else: f = np.array([pdf.xfxQ2(4, x, muf2) / x for x in out["interpolation_xgrid"]]) -# check yadism SV is proportional splitting functions -# print("check yadism SV") -# if mode_g: -# print((pqgg @ f)/(loq @ pqg.T @ f)) -# else: -# print((pqqc @ f)/(loc @ pqq.T @ f)) + +def check_yad_sv(plot=False): + """Check if yadism SV is proportional to splitting functions.""" + print("check yadism SV") + # collect data + ops = { + "nlosv": None, + "nnlosv1": None, + "nnlosv2": None, + } + if mode_g: + ops["nlosv"] = (nlosvg, loq @ pqg0) + ops["nnlosv1"] = ( + nnlosv1g, + loq @ pqg1 + nlog @ (pgg0 - beta0I) + nloq @ pqg0, + ) + ops["nnlosv2"] = ( + nnlosv2g, + 1.0 / 2.0 * loq @ pqq0 @ pqg0 + 1.0 / 2.0 * (loq @ pqg0 @ (pgg0 - beta0I)), + ) + else: + ops["nlosv"] = (nlosvc, loc @ pqq0) + ops["nnlosv1"] = ( + nnlosv1c, + loc @ pqq1 + nloc @ (pqq0 - beta0I) + nlog @ pgq0, + ) + ops["nnlosv2"] = ( + nnlosv2c, + 1.0 / 2.0 * loq @ pqg0 @ pgq0 + 1.0 / 2.0 * (loc @ pqq0 @ (pqq0 - beta0I)), + ) + # print check + for k, v in ops.items(): + if v is None: + continue + print("ratio", k) + print((v[0] @ f) / (v[1] @ f)) + # show actual curves? + if plot: + plt.title(f"{k}@{pdf.set().name}") + plt.plot((v[0] @ f), label="yad") + plt.plot((v[1] @ f), label="ana") + plt.legend() + plt.show() + + +check_yad_sv(True) # extract a by comparing C against yadism -# print("check C result") -# if mode_g: -# print(cc1/((ac*(nlog + pqgg*L)) @ f)) -# else: -# print(cc1/((loc + ac*(nloc + pqqc*L)) @ f)) +def check_c_res(): + print("check C result") + if mode_g: + op = ac * (nlog + nlosvg * L) + else: + op = loc + ac * (nloc + nlosvc * L) + print(cc1 / (op @ f)) + + +# check_c_res() # check B EKO is K # print("check B EKO result") @@ -181,26 +297,25 @@ # plt.show() # check difference between B and C -print("check B - C") -diff_grid = bb1 - cc1 -# plt.plot(bb1, label="b-grid") -# plt.plot(cc1, label="c-grid") -# nloc *= 1.9 -# nlog *= 1.14 -if mode_g: - # plt.plot(a**2 * L * ((nloc @ pqg.T @ f)), label="pqg") - # plt.plot(a**2 * L * ((nlog @ pgg.T @ f)), label="pgg") - diff_ana = (nlog * (ab - ac) + ac * ab * L * ((nloq @ pqg.T) + (nlog @ pgg.T))) @ f -else: - # plt.plot(a**2 * L * ((nloc @ pqq.T @ f)), label="pqq") - # plt.plot(a**2 * L * ((nlog @ pgq.T @ f)), label="pgq") - diff_ana = (nloc * (ab - ac) + ab * ac * L * ((nloc @ pqq.T) + (nlog @ pgq.T))) @ f -# print(diff_grid/cc1) -print(bb1 / cc1) -print(diff_grid / diff_ana) -plt.plot(diff_grid, label="grid") -plt.plot(diff_ana, label="ana") -# plt.plot((diff_grid/bb1-1)*100, label="rel err. diff to b") -# plt.plot(diff_grid/diff_ana, label="ratio") -plt.legend() -plt.show() +# print("check B - C") +# if order_mode == Order.nlo or order_mode == Order.nnlo: +# diff_grid = bb1 - cc1 +# # plt.plot(bb1, label="b-grid") +# # plt.plot(cc1, label="c-grid") +# if mode_g: +# # plt.plot(a**2 * L * ((nloc @ pqg0 @ f)), label="pqg0") +# # plt.plot(a**2 * L * ((nlog @ pgg0 @ f)), label="pgg0") +# diff_ana = (nlog * (ab - ac) + ac * ab * L * ((nloq @ pqg0) + (nlog @ pgg0))) @ f +# else: +# # plt.plot(a**2 * L * ((nloc @ pqq0 @ f)), label="pqq0") +# # plt.plot(a**2 * L * ((nlog @ pgq0 @ f)), label="pgq0") +# diff_ana = (nloc * (ab - ac) + ab * ac * L * ((nloc @ pqq0) + (nlog @ pgq0))) @ f +# # print(diff_grid/cc1) +# print(bb1 / cc1) +# print(diff_grid / diff_ana) +# # plt.plot(diff_grid, label="grid") +# # plt.plot(diff_ana, label="ana") +# # plt.plot((diff_grid/bb1-1)*100, label="rel err. diff to b") +# # plt.plot(diff_grid/diff_ana, label="ratio") +# # plt.legend() +# # plt.show() From 74dcaea4b0448c576e2ab4c4aed1f2e81c731b23 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Wed, 2 Nov 2022 10:43:04 +0100 Subject: [PATCH 09/24] Add things to be compatible with eko --- debug3.py | 58 +++++++++++++++++++++++++++++-------------------------- 1 file changed, 31 insertions(+), 27 deletions(-) diff --git a/debug3.py b/debug3.py index f08377f0..d3d1ee4a 100644 --- a/debug3.py +++ b/debug3.py @@ -1,18 +1,18 @@ -# -*- coding: utf-8 -*- -from enum import Enum import sys +from enum import Enum + +import eko import lhapdf -import numpy as np import matplotlib.pyplot as plt +import numpy as np import pineappl import yadism import yaml -import eko -from eko import interpolation -from eko import output -from eko import strong_coupling as sc +from eko import compatibility +from eko import couplings as sc +from eko import interpolation, output +from eko.beta import beta_qcd_as2 # find the substitute from ekomark.apply import apply_pdf -from eko.beta import beta_0 # from eko.output.struct import EKO from yadism.coefficient_functions.splitting_functions import lo, nlo @@ -36,7 +36,7 @@ class Order(Enum): order_mode = Order.nonlo if len(sys.argv) > 2: - om = sys.argv[1].strip().lower() + om = sys.argv[2].strip().lower() if om == "nlo": order_mode = Order.nlo elif om == "nnlo": @@ -47,19 +47,21 @@ class Order(Enum): pdf = lhapdf.mkPDF("gonly", 0) else: pdf = lhapdf.mkPDF("conly", 0) - # load FK tables if order_mode == Order.nonlo: - b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/no-nlo-test.pineappl.lz4") - c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/no-nlo-test.pineappl.lz4") + b1 = pineappl.fk_table.FkTable.read( + "../data/fktables/2205/no-nlo-test.pineappl.lz4" + ) + c1 = pineappl.fk_table.FkTable.read( + "../data/fktables/3205/no-nlo-test.pineappl.lz4" + ) elif order_mode == Order.nlo: - b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/test.pineappl.lz4") - c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/test.pineappl.lz4") + b1 = pineappl.fk_table.FkTable.read("../data/fktables/2205/test.pineappl.lz4") + c1 = pineappl.fk_table.FkTable.read("../data/fktables/3205/test.pineappl.lz4") elif order_mode == Order.nnlo: - b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/test-nnlo.pineappl.lz4") - c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/test-nnlo.pineappl.lz4") -out = yadism.output.Output.load_tar("test0.tar") - + b1 = pineappl.fk_table.FkTable.read("../data/fktables/2205/test-nnlo.pineappl.lz4") + c1 = pineappl.fk_table.FkTable.read("../data/fktables/3205/test-nnlo.pineappl.lz4") +out = yadism.output.Output.load_tar("../test0.tar") # compute predictions from FK bb1 = b1.convolute_with_one(2212, pdf.xfxQ2) cc1 = c1.convolute_with_one(2212, pdf.xfxQ2) @@ -156,23 +158,25 @@ class Order(Enum): ) # compute splitting functions -interp = interpolation.InterpolatorDispatcher.from_dict(out, mode_N=False) +mygrid = eko.interpolation.XGrid(out["interpolation_xgrid"]) +interp = interpolation.InterpolatorDispatcher(mygrid, 4) pqq0 = convolute_operator(lo.pqq(5), interp)[0].T pqg0 = convolute_operator(lo.pqg(5), interp)[0].T / 10 # * 55 / 100 pgq0 = convolute_operator(nlo.pgq0(5), interp)[0].T pgg0 = convolute_operator(nlo.pgg0(5), interp)[0].T pqq1 = convolute_operator(nlo.pqq1(5), interp)[0].T pqg1 = convolute_operator(nlo.pqg1(5), interp)[0].T -beta0I = beta_0(5) * np.eye(pqg0.shape[0]) +beta0I = beta_qcd_as2(5) * np.eye(pqg0.shape[0]) q2 = 300.0 muf2 = q2 * 2.0**2 # load theory card with open("data/theory_cards/2205.yaml") as fd: - tc = yaml.safe_load(fd) -astrong = sc.StrongCoupling.from_dict(tc) + old_tc = yaml.safe_load(fd) +tc = eko.compatibility.update_theory(old_tc) +astrong = sc.Couplings.from_dict(tc) # actually, there are two alpha_s values in the game! -ac = astrong.a_s(muf2) +# ac = astrong.a_s(muf2) ab = astrong.a_s(q2) L = np.log(1.0 / 2.0**2) @@ -192,10 +196,10 @@ class Order(Enum): nnloq = 2 * (3 * nnlos + 2 * nnloc) # construct K -Kqq = np.eye(len(out["interpolation_xgrid"])) + ac * L * pqq0 -Kgg = np.eye(len(out["interpolation_xgrid"])) + ac * L * pgg0 -Kqg = ac * L * pqg0 -Kgq = ac * L * pgq0 +Kqq = np.eye(len(out["interpolation_xgrid"])) + ab * L * pqq0 +Kgg = np.eye(len(out["interpolation_xgrid"])) + ab * L * pgg0 +Kqg = ab * L * pqg0 +Kgq = ab * L * pgq0 # compute PDFs if mode_g: From 972315fbf9799319d04a656550889267b81f997c Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Wed, 2 Nov 2022 12:16:42 +0100 Subject: [PATCH 10/24] Fix interpolator --- debug3.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/debug3.py b/debug3.py index d3d1ee4a..a24d4d3d 100644 --- a/debug3.py +++ b/debug3.py @@ -159,7 +159,9 @@ class Order(Enum): # compute splitting functions mygrid = eko.interpolation.XGrid(out["interpolation_xgrid"]) -interp = interpolation.InterpolatorDispatcher(mygrid, 4) +interp = interpolation.InterpolatorDispatcher( + mygrid, out["interpolation_polynomial_degree"], mode_N=False +) pqq0 = convolute_operator(lo.pqq(5), interp)[0].T pqg0 = convolute_operator(lo.pqg(5), interp)[0].T / 10 # * 55 / 100 pgq0 = convolute_operator(nlo.pgq0(5), interp)[0].T @@ -252,19 +254,19 @@ def check_yad_sv(plot=False): plt.show() -check_yad_sv(True) +# check_yad_sv(True) # extract a by comparing C against yadism def check_c_res(): print("check C result") if mode_g: - op = ac * (nlog + nlosvg * L) + op = ab * (nlog + nlosvg * L) else: - op = loc + ac * (nloc + nlosvc * L) + op = loc + ab * (nloc + nlosvc * L) print(cc1 / (op @ f)) -# check_c_res() +check_c_res() # check B EKO is K # print("check B EKO result") From 9815399907b30dbb63d12f0b261c4201516c8eb3 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Wed, 2 Nov 2022 15:25:16 +0100 Subject: [PATCH 11/24] Completed compatibility with new eko --- debug3.py | 92 ++++++++++++++++++++++++++++--------------------------- 1 file changed, 47 insertions(+), 45 deletions(-) diff --git a/debug3.py b/debug3.py index a24d4d3d..a7c5d4fd 100644 --- a/debug3.py +++ b/debug3.py @@ -12,12 +12,15 @@ from eko import couplings as sc from eko import interpolation, output from eko.beta import beta_qcd_as2 # find the substitute +from eko.output import legacy from ekomark.apply import apply_pdf # from eko.output.struct import EKO from yadism.coefficient_functions.splitting_functions import lo, nlo from yadism.esf.conv import convolute_operator +from pineko import ekompatibility + # mode can be gluon only or charm only mode_g = True if len(sys.argv) > 1: @@ -267,58 +270,57 @@ def check_c_res(): check_c_res() - # check B EKO is K -# print("check B EKO result") -# ekob = output.Output.load_tar("data/ekos/2205/test.tar") -# fb = apply_pdf(ekob, pdf) -# if mode_g: -# gb = Kgg @ f -# cb = Kqg @ f -# else: -# gb = Kgq @ f -# cb = Kqq @ f -# print(fb[300.0]["pdfs"][21] / gb) -# print(fb[300.0]["pdfs"][4] / cb) +print("check B EKO result") +ekob = legacy.load_tar("../data/ekos/2205/test.tar") +fb = apply_pdf(ekob, pdf) +if mode_g: + gb = Kgg @ f + cb = Kqg @ f +else: + gb = Kgq @ f + cb = Kqq @ f +print(fb[300.0]["pdfs"][21] / gb) +print(fb[300.0]["pdfs"][4] / cb) # check C EKO is identity -# print("check C EKO result") -# ekoc = output.Output.load_tar("data/ekos/3205/test.tar") -# fc = apply_pdf(ekoc, pdf) -# if mode_g: -# print(fc[2.**2 * 300.]["pdfs"][21]/f) -# print(fc[2.**2 * 300.]["pdfs"][4] - np.zeros_like(f)) -# else: -# print(fc[2.**2 * 300.]["pdfs"][4]/f) -# print(fc[2.**2 * 300.]["pdfs"][21] - np.zeros_like(f)) +print("check C EKO result") +ekoc = legacy.load_tar("../data/ekos/3205/test.tar") +fc = apply_pdf(ekoc, pdf) +if mode_g: + print(fc[2.0**2 * 300.0]["pdfs"][21] / f) + print(fc[2.0**2 * 300.0]["pdfs"][4] - np.zeros_like(f)) +else: + print(fc[2.0**2 * 300.0]["pdfs"][4] / f) + print(fc[2.0**2 * 300.0]["pdfs"][21] - np.zeros_like(f)) # check B result -# print("check B result") -# if mode_g: -# print(bb1 / (((loq + ab * nloq) @ Kqg + (ab * nlog) @ Kgg) @ f)) -# plt.plot(bb1) -# plt.plot((((loq + ab * nloq) @ Kqg + (ab * nlog) @ Kgg) @ f)) -# else: -# print(bb1 / (((loc + ab * nloc) @ Kqq + (ab * nlog) @ Kgq) @ f)) -# plt.show() +print("check B result") +if mode_g: + print(bb1 / (((loq + ab * nloq) @ Kqg + (ab * nlog) @ Kgg) @ f)) + plt.plot(bb1) + plt.plot(((loq + ab * nloq) @ Kqg + (ab * nlog) @ Kgg) @ f) +else: + print(bb1 / (((loc + ab * nloc) @ Kqq + (ab * nlog) @ Kgq) @ f)) +plt.show() # check difference between B and C -# print("check B - C") -# if order_mode == Order.nlo or order_mode == Order.nnlo: -# diff_grid = bb1 - cc1 -# # plt.plot(bb1, label="b-grid") -# # plt.plot(cc1, label="c-grid") -# if mode_g: -# # plt.plot(a**2 * L * ((nloc @ pqg0 @ f)), label="pqg0") -# # plt.plot(a**2 * L * ((nlog @ pgg0 @ f)), label="pgg0") -# diff_ana = (nlog * (ab - ac) + ac * ab * L * ((nloq @ pqg0) + (nlog @ pgg0))) @ f -# else: -# # plt.plot(a**2 * L * ((nloc @ pqq0 @ f)), label="pqq0") -# # plt.plot(a**2 * L * ((nlog @ pgq0 @ f)), label="pgq0") -# diff_ana = (nloc * (ab - ac) + ab * ac * L * ((nloc @ pqq0) + (nlog @ pgq0))) @ f -# # print(diff_grid/cc1) -# print(bb1 / cc1) -# print(diff_grid / diff_ana) +print("check B - C") +if order_mode == Order.nlo or order_mode == Order.nnlo: + diff_grid = bb1 - cc1 + plt.plot(bb1, label="b-grid") + plt.plot(cc1, label="c-grid") + if mode_g: + # plt.plot(ab**2 * L * ((nloc @ pqg0 @ f)), label="pqg0") + # plt.plot(ab**2 * L * ((nlog @ pgg0 @ f)), label="pgg0") + diff_ana = (ab * ab * L * ((nloq @ pqg0) + (nlog @ pgg0))) @ f + else: + # plt.plot(a**2 * L * ((nloc @ pqq0 @ f)), label="pqq0") + # plt.plot(a**2 * L * ((nlog @ pgq0 @ f)), label="pgq0") + diff_ana = (ab * ab * L * ((nloc @ pqq0) + (nlog @ pgq0))) @ f + print(diff_grid / cc1) + print(bb1 / cc1) + print(diff_grid / diff_ana) # # plt.plot(diff_grid, label="grid") # # plt.plot(diff_ana, label="ana") # # plt.plot((diff_grid/bb1-1)*100, label="rel err. diff to b") From 037ed1c74941547c1cb07e12a489ed7ebac5437e Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Fri, 4 Nov 2022 11:21:13 +0100 Subject: [PATCH 12/24] Improve debug3 --- debug3.py | 30 +++++++++++++++++++++--------- 1 file changed, 21 insertions(+), 9 deletions(-) diff --git a/debug3.py b/debug3.py index a7c5d4fd..d22edca6 100644 --- a/debug3.py +++ b/debug3.py @@ -205,6 +205,14 @@ class Order(Enum): Kgg = np.eye(len(out["interpolation_xgrid"])) + ab * L * pgg0 Kqg = ab * L * pqg0 Kgq = ab * L * pgq0 +# We believe to have a minus sign problem so: +minus_sign = True +if minus_sign: + Kqg = -Kqg + Kgq = -Kgq + Kqq = np.eye(len(out["interpolation_xgrid"])) - ab * L * pqq0 + Kgg = np.eye(len(out["interpolation_xgrid"])) - ab * L * pgg0 + # compute PDFs if mode_g: @@ -298,18 +306,22 @@ def check_c_res(): print("check B result") if mode_g: print(bb1 / (((loq + ab * nloq) @ Kqg + (ab * nlog) @ Kgg) @ f)) - plt.plot(bb1) - plt.plot(((loq + ab * nloq) @ Kqg + (ab * nlog) @ Kgg) @ f) + plt.plot(bb1, label="b-fk") + plt.plot(((loq + ab * nloq) @ (Kqg) + (ab * nlog) @ (Kgg)) @ f, label="ana_pred") else: print(bb1 / (((loc + ab * nloc) @ Kqq + (ab * nlog) @ Kgq) @ f)) + plt.plot(bb1, label="b-fk") + plt.plot(((loc + ab * nloc) @ Kqq + (ab * nlog) @ Kgq) @ f, label="ana_pred") +plt.legend() plt.show() - # check difference between B and C print("check B - C") if order_mode == Order.nlo or order_mode == Order.nnlo: diff_grid = bb1 - cc1 plt.plot(bb1, label="b-grid") plt.plot(cc1, label="c-grid") + plt.legend() + plt.show() if mode_g: # plt.plot(ab**2 * L * ((nloc @ pqg0 @ f)), label="pqg0") # plt.plot(ab**2 * L * ((nlog @ pgg0 @ f)), label="pgg0") @@ -321,9 +333,9 @@ def check_c_res(): print(diff_grid / cc1) print(bb1 / cc1) print(diff_grid / diff_ana) -# # plt.plot(diff_grid, label="grid") -# # plt.plot(diff_ana, label="ana") -# # plt.plot((diff_grid/bb1-1)*100, label="rel err. diff to b") -# # plt.plot(diff_grid/diff_ana, label="ratio") -# # plt.legend() -# # plt.show() + plt.plot(diff_grid, label="grid") + plt.plot(diff_ana, label="ana") + # # plt.plot((diff_grid/bb1-1)*100, label="rel err. diff to b") + # # plt.plot(diff_grid/diff_ana, label="ratio") + plt.legend() + plt.show() From 4c81c88b6cceda49e06715dd7ff5eb5cd0a1531e Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Mon, 7 Nov 2022 10:54:30 +0100 Subject: [PATCH 13/24] Complete test on schemes --- debug3.py | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/debug3.py b/debug3.py index d22edca6..e2ae1e98 100644 --- a/debug3.py +++ b/debug3.py @@ -45,6 +45,12 @@ class Order(Enum): elif om == "nnlo": order_mode = Order.nnlo +signed = False +if len(sys.argv) > 3: + sig = sys.argv[3].strip().lower() + if sig == "signed": + signed = True + # load respective PDF set if mode_g: pdf = lhapdf.mkPDF("gonly", 0) @@ -206,8 +212,7 @@ class Order(Enum): Kqg = ab * L * pqg0 Kgq = ab * L * pgq0 # We believe to have a minus sign problem so: -minus_sign = True -if minus_sign: +if signed: Kqg = -Kqg Kgq = -Kgq Kqq = np.eye(len(out["interpolation_xgrid"])) - ab * L * pqq0 @@ -326,16 +331,25 @@ def check_c_res(): # plt.plot(ab**2 * L * ((nloc @ pqg0 @ f)), label="pqg0") # plt.plot(ab**2 * L * ((nlog @ pgg0 @ f)), label="pgg0") diff_ana = (ab * ab * L * ((nloq @ pqg0) + (nlog @ pgg0))) @ f + if signed: + diff_ana = ( + -2 * ab * L * (loq @ pqg0) + - ab * ab * L * ((nloq @ pqg0) + (nlog @ pgg0)) + ) @ f else: # plt.plot(a**2 * L * ((nloc @ pqq0 @ f)), label="pqq0") # plt.plot(a**2 * L * ((nlog @ pgq0 @ f)), label="pgq0") diff_ana = (ab * ab * L * ((nloc @ pqq0) + (nlog @ pgq0))) @ f - print(diff_grid / cc1) - print(bb1 / cc1) - print(diff_grid / diff_ana) + if signed: + diff_ana = ( + -2 * ab * L * (loc @ pqq0) + - ab * ab * L * ((nloc @ pqq0) + (nlog @ pgq0)) + ) @ f + print(((diff_grid / diff_ana) - 1.0) * 100) plt.plot(diff_grid, label="grid") plt.plot(diff_ana, label="ana") - # # plt.plot((diff_grid/bb1-1)*100, label="rel err. diff to b") - # # plt.plot(diff_grid/diff_ana, label="ratio") + plt.legend() + plt.show() + plt.plot(((diff_grid[1:-5] / diff_ana[1:-5]) - 1.0) * 100, label="ratio") plt.legend() plt.show() From 9f6a1b0132c1830260f8c795ca213ffbbef20743 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Wed, 2 Nov 2022 10:48:08 +0100 Subject: [PATCH 14/24] Fix iterate call in opcards --- src/pineko/theory.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/pineko/theory.py b/src/pineko/theory.py index 05a236e7..5242aa50 100644 --- a/src/pineko/theory.py +++ b/src/pineko/theory.py @@ -239,7 +239,9 @@ def opcards(self): """Write operator cards.""" tcard = theory_card.load(self.theory_id) self.operator_cards_path.mkdir(exist_ok=True) - self.iterate(self.opcard, xif=tcard["XIF"]) + self.iterate( + self.opcard, xif=tcard["XIF"], tcard_path=theory_card.path(self.theory_id) + ) def load_operator_card(self, name): """Read current operator card. From 8a030a6df07f4d60cdb66dc3c9edb196be81411d Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Wed, 2 Nov 2022 12:20:35 +0100 Subject: [PATCH 15/24] Fix dump of eko --- src/pineko/theory.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pineko/theory.py b/src/pineko/theory.py index 5242aa50..cbb4ea31 100644 --- a/src/pineko/theory.py +++ b/src/pineko/theory.py @@ -320,7 +320,7 @@ def eko(self, name, _grid, tcard): logger.info("Start computation of %s", name) start_time = time.perf_counter() ops = eko.run_dglap(theory_card=tcard, operators_card=ocard) - ops.dump_tar(eko_filename) + ops.deepcopy(eko_filename) logger.info( "Finished computation of %s - took %f s", name, From 15177444bfef6655993e1d0b7bc86c6fd1472f42 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Wed, 2 Nov 2022 12:31:33 +0100 Subject: [PATCH 16/24] Fix eko dump again --- src/pineko/theory.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pineko/theory.py b/src/pineko/theory.py index cbb4ea31..7576460f 100644 --- a/src/pineko/theory.py +++ b/src/pineko/theory.py @@ -320,7 +320,7 @@ def eko(self, name, _grid, tcard): logger.info("Start computation of %s", name) start_time = time.perf_counter() ops = eko.run_dglap(theory_card=tcard, operators_card=ocard) - ops.deepcopy(eko_filename) + eko.output.legacy.dump_tar(ops, eko_filename) logger.info( "Finished computation of %s - took %f s", name, From 9d5c260b37f2c52fb1e4357740cf9bad5d9a9388 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Wed, 2 Nov 2022 13:49:00 +0100 Subject: [PATCH 17/24] Remove inputpids and outputpids from op card --- benchmarks/bench_evolve.py | 12 +++++++++--- src/pineko/evolve.py | 10 +++++++++- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/benchmarks/bench_evolve.py b/benchmarks/bench_evolve.py index b43a5fba..018d223e 100644 --- a/benchmarks/bench_evolve.py +++ b/benchmarks/bench_evolve.py @@ -32,9 +32,15 @@ def benchmark_write_operator_card_from_file(tmp_path, test_files, test_configs): assert np.allclose(myopcard["rotations"]["xgrid"], x_grid) assert np.allclose(myopcard["rotations"]["targetgrid"], x_grid) assert np.allclose(myopcard["rotations"]["inputgrid"], x_grid) - assert np.allclose(myopcard["rotations"]["inputpids"], pineko.evolve.DEFAULT_PIDS) - assert np.allclose(myopcard["rotations"]["targetpids"], pineko.evolve.DEFAULT_PIDS) - + # I am now testing if the keys that I removed are missing as they should + with pytest.raises(KeyError): + myopcard["rotations"]["inputpids"] + with pytest.raises(KeyError): + myopcard["rotations"]["targetpids"] + with pytest.raises(KeyError): + myopcard["inputpids"] + with pytest.raises(KeyError): + myopcard["targetpids"] wrong_pine_path = test_files / "data/grids/208/HERA_CC_318GEV_EM_wrong.pineappl.lz4" with pytest.raises(FileNotFoundError): _ = pineko.evolve.write_operator_card_from_file( diff --git a/src/pineko/evolve.py b/src/pineko/evolve.py index ea81cb66..ed190a5c 100644 --- a/src/pineko/evolve.py +++ b/src/pineko/evolve.py @@ -108,8 +108,16 @@ def provide_if_missing(key, default, card=operators_card): provide_if_missing("targetgrid", operators_card["interpolation_xgrid"]) provide_if_missing("inputpids", DEFAULT_PIDS) provide_if_missing("targetpids", DEFAULT_PIDS) - _, new_operators_card = eko.compatibility.update(tcard, operators_card) + + def remove_useless_entries(key, card=new_operators_card): + if key in card: + card.pop(key) + if key in card["rotations"]: + card["rotations"].pop(key) + + remove_useless_entries("inputpids") + remove_useless_entries("targetpids") with open(card_path, "w", encoding="UTF-8") as f: yaml.safe_dump(new_operators_card, f) return x_grid, q2_grid From f01e7d7db5dfe0f770eb60c79240ac235ef7f168 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Wed, 2 Nov 2022 15:44:46 +0100 Subject: [PATCH 18/24] Revert "Remove inputpids and outputpids from op card" This reverts commit 359c4547919cd6109611e816e4b8718403c4fb01. --- benchmarks/bench_evolve.py | 12 +++--------- src/pineko/evolve.py | 10 +--------- 2 files changed, 4 insertions(+), 18 deletions(-) diff --git a/benchmarks/bench_evolve.py b/benchmarks/bench_evolve.py index 018d223e..b43a5fba 100644 --- a/benchmarks/bench_evolve.py +++ b/benchmarks/bench_evolve.py @@ -32,15 +32,9 @@ def benchmark_write_operator_card_from_file(tmp_path, test_files, test_configs): assert np.allclose(myopcard["rotations"]["xgrid"], x_grid) assert np.allclose(myopcard["rotations"]["targetgrid"], x_grid) assert np.allclose(myopcard["rotations"]["inputgrid"], x_grid) - # I am now testing if the keys that I removed are missing as they should - with pytest.raises(KeyError): - myopcard["rotations"]["inputpids"] - with pytest.raises(KeyError): - myopcard["rotations"]["targetpids"] - with pytest.raises(KeyError): - myopcard["inputpids"] - with pytest.raises(KeyError): - myopcard["targetpids"] + assert np.allclose(myopcard["rotations"]["inputpids"], pineko.evolve.DEFAULT_PIDS) + assert np.allclose(myopcard["rotations"]["targetpids"], pineko.evolve.DEFAULT_PIDS) + wrong_pine_path = test_files / "data/grids/208/HERA_CC_318GEV_EM_wrong.pineappl.lz4" with pytest.raises(FileNotFoundError): _ = pineko.evolve.write_operator_card_from_file( diff --git a/src/pineko/evolve.py b/src/pineko/evolve.py index ed190a5c..ea81cb66 100644 --- a/src/pineko/evolve.py +++ b/src/pineko/evolve.py @@ -108,16 +108,8 @@ def provide_if_missing(key, default, card=operators_card): provide_if_missing("targetgrid", operators_card["interpolation_xgrid"]) provide_if_missing("inputpids", DEFAULT_PIDS) provide_if_missing("targetpids", DEFAULT_PIDS) - _, new_operators_card = eko.compatibility.update(tcard, operators_card) - - def remove_useless_entries(key, card=new_operators_card): - if key in card: - card.pop(key) - if key in card["rotations"]: - card["rotations"].pop(key) - remove_useless_entries("inputpids") - remove_useless_entries("targetpids") + _, new_operators_card = eko.compatibility.update(tcard, operators_card) with open(card_path, "w", encoding="UTF-8") as f: yaml.safe_dump(new_operators_card, f) return x_grid, q2_grid From 8a74ca8f4260e9c7b0a89d0c603e93eade796878 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Wed, 2 Nov 2022 15:45:03 +0100 Subject: [PATCH 19/24] Revert "Fix eko dump again" This reverts commit 34bfc72eb4ae943b3e27fc9172128fe0d7730192. --- src/pineko/theory.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pineko/theory.py b/src/pineko/theory.py index 7576460f..cbb4ea31 100644 --- a/src/pineko/theory.py +++ b/src/pineko/theory.py @@ -320,7 +320,7 @@ def eko(self, name, _grid, tcard): logger.info("Start computation of %s", name) start_time = time.perf_counter() ops = eko.run_dglap(theory_card=tcard, operators_card=ocard) - eko.output.legacy.dump_tar(ops, eko_filename) + ops.deepcopy(eko_filename) logger.info( "Finished computation of %s - took %f s", name, From 8367f9919604eb84e2f3b424679040e22cff9ab3 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Wed, 2 Nov 2022 15:45:15 +0100 Subject: [PATCH 20/24] Revert "Fix dump of eko" This reverts commit 3effc686c0b446d508734b23d25865dd5c3c6325. --- src/pineko/theory.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pineko/theory.py b/src/pineko/theory.py index cbb4ea31..5242aa50 100644 --- a/src/pineko/theory.py +++ b/src/pineko/theory.py @@ -320,7 +320,7 @@ def eko(self, name, _grid, tcard): logger.info("Start computation of %s", name) start_time = time.perf_counter() ops = eko.run_dglap(theory_card=tcard, operators_card=ocard) - ops.deepcopy(eko_filename) + ops.dump_tar(eko_filename) logger.info( "Finished computation of %s - took %f s", name, From 2820a2ac989b3719069ec64cda83f33f21803d7a Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Thu, 10 Nov 2022 12:31:06 +0100 Subject: [PATCH 21/24] Allow nnlo checks --- debug3.py | 189 ++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 148 insertions(+), 41 deletions(-) diff --git a/debug3.py b/debug3.py index e2ae1e98..31beed2d 100644 --- a/debug3.py +++ b/debug3.py @@ -11,7 +11,7 @@ from eko import compatibility from eko import couplings as sc from eko import interpolation, output -from eko.beta import beta_qcd_as2 # find the substitute +from eko.beta import beta_qcd_as2 from eko.output import legacy from ekomark.apply import apply_pdf @@ -58,18 +58,17 @@ class Order(Enum): pdf = lhapdf.mkPDF("conly", 0) # load FK tables if order_mode == Order.nonlo: - b1 = pineappl.fk_table.FkTable.read( - "../data/fktables/2205/no-nlo-test.pineappl.lz4" - ) - c1 = pineappl.fk_table.FkTable.read( - "../data/fktables/3205/no-nlo-test.pineappl.lz4" - ) + b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/no-nlo-test.pineappl.lz4") + c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/no-nlo-test.pineappl.lz4") elif order_mode == Order.nlo: - b1 = pineappl.fk_table.FkTable.read("../data/fktables/2205/test.pineappl.lz4") - c1 = pineappl.fk_table.FkTable.read("../data/fktables/3205/test.pineappl.lz4") + b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/test.pineappl.lz4") + c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/test.pineappl.lz4") elif order_mode == Order.nnlo: - b1 = pineappl.fk_table.FkTable.read("../data/fktables/2205/test-nnlo.pineappl.lz4") - c1 = pineappl.fk_table.FkTable.read("../data/fktables/3205/test-nnlo.pineappl.lz4") + # b1 = pineappl.fk_table.FkTable.read("data/fktables/2405/test-nnlo.pineappl.lz4") #this is the one with the wrong sign + c1 = pineappl.fk_table.FkTable.read("data/fktables/3405/test-nnlo.pineappl.lz4") + b1 = pineappl.fk_table.FkTable.read( + "data/fktables/2505/test-nnlo.pineappl.lz4" + ) # this is the one produced with the tentative solution of the sign out = yadism.output.Output.load_tar("../test0.tar") # compute predictions from FK bb1 = b1.convolute_with_one(2212, pdf.xfxQ2) @@ -176,14 +175,22 @@ class Order(Enum): pgq0 = convolute_operator(nlo.pgq0(5), interp)[0].T pgg0 = convolute_operator(nlo.pgg0(5), interp)[0].T pqq1 = convolute_operator(nlo.pqq1(5), interp)[0].T +# pgg1 = convolute_operator(nlo.pgg1(5), interp)[0].T +pgg1 = pgg0 pqg1 = convolute_operator(nlo.pqg1(5), interp)[0].T +# pgq1 = convolute_operator(nlo.pgq1(5), interp)[0].T +pgq1 = pgq0 beta0I = beta_qcd_as2(5) * np.eye(pqg0.shape[0]) - +beta0 = beta_qcd_as2(5) q2 = 300.0 muf2 = q2 * 2.0**2 # load theory card -with open("data/theory_cards/2205.yaml") as fd: - old_tc = yaml.safe_load(fd) +if order_mode == Order.nnlo: + with open("data/theory_cards/2405.yaml") as fd: + old_tc = yaml.safe_load(fd) +else: + with open("data/theory_cards/2205.yaml") as fd: + old_tc = yaml.safe_load(fd) tc = eko.compatibility.update_theory(old_tc) astrong = sc.Couplings.from_dict(tc) # actually, there are two alpha_s values in the game! @@ -224,8 +231,7 @@ class Order(Enum): f = np.array([pdf.xfxQ2(21, x, muf2) / x for x in out["interpolation_xgrid"]]) else: f = np.array([pdf.xfxQ2(4, x, muf2) / x for x in out["interpolation_xgrid"]]) - - +# Todo: Correct this part at nnlo def check_yad_sv(plot=False): """Check if yadism SV is proportional to splitting functions.""" print("check yadism SV") @@ -239,21 +245,24 @@ def check_yad_sv(plot=False): ops["nlosv"] = (nlosvg, loq @ pqg0) ops["nnlosv1"] = ( nnlosv1g, - loq @ pqg1 + nlog @ (pgg0 - beta0I) + nloq @ pqg0, + loq @ pqg1 + nlog @ pgg0 + nloq @ pqg0, ) ops["nnlosv2"] = ( nnlosv2g, - 1.0 / 2.0 * loq @ pqq0 @ pqg0 + 1.0 / 2.0 * (loq @ pqg0 @ (pgg0 - beta0I)), + (1.0 / 2.0) * loq @ pqq0 @ pqg0 + + (-beta0) * (1.0 / 2.0) * loq @ pqq0 + + (1.0 / 2.0) * (loq @ pqg0 @ pgg0) + + (-beta0) * (1.0 / 2.0) * (loq @ pqg0), ) else: ops["nlosv"] = (nlosvc, loc @ pqq0) ops["nnlosv1"] = ( nnlosv1c, - loc @ pqq1 + nloc @ (pqq0 - beta0I) + nlog @ pgq0, + loc @ pqq1 + nloc @ pqq0 + nlog @ pgq0, ) ops["nnlosv2"] = ( nnlosv2c, - 1.0 / 2.0 * loq @ pqg0 @ pgq0 + 1.0 / 2.0 * (loc @ pqq0 @ (pqq0 - beta0I)), + 1.0 / 2.0 * loq @ pqg0 @ pgq0 + 1.0 / 2.0 * (loc @ pqq0 @ (pqq0 + beta0I)), ) # print check for k, v in ops.items(): @@ -276,16 +285,30 @@ def check_yad_sv(plot=False): def check_c_res(): print("check C result") if mode_g: - op = ab * (nlog + nlosvg * L) + if order_mode == Order.nlo or order_mode == Order.nonlo: + op = ab * (nlog + nlosvg * L) + if order_mode == Order.nnlo: + op_nlo = ab * (nlog + nlosvg * L) + op_nnlo = ab * ab * (nnlog + L * nnlosv1g + L * L * nnlosv2g) + op = op_nlo + op_nnlo else: - op = loc + ab * (nloc + nlosvc * L) + if order_mode == Order.nlo or order_mode == Order.nonlo: + op = loc + ab * (nloc + nlosvc * L) + if order_mode == Order.nnlo: + op_nlo = loc + ab * (nloc + nlosvc * L) + op_nnlo = ab * ab * (nnlog + L * nnlosv1c + L * L * nnlosv2c) + op = op_nlo + op_nnlo + print(cc1 / (op @ f)) check_c_res() # check B EKO is K print("check B EKO result") -ekob = legacy.load_tar("../data/ekos/2205/test.tar") +if order_mode == Order.nnlo: + ekob = output.EKO.load("data/ekos/2405/test-nnlo.tar") +else: + ekob = legacy.load_tar("data/ekos/2205/test.tar") fb = apply_pdf(ekob, pdf) if mode_g: gb = Kgg @ f @@ -298,7 +321,10 @@ def check_c_res(): # check C EKO is identity print("check C EKO result") -ekoc = legacy.load_tar("../data/ekos/3205/test.tar") +if order_mode == Order.nnlo: + ekoc = output.EKO.load("data/ekos/3405/test-nnlo.tar") +else: + ekoc = legacy.load_tar("data/ekos/3205/test.tar") fc = apply_pdf(ekoc, pdf) if mode_g: print(fc[2.0**2 * 300.0]["pdfs"][21] / f) @@ -311,14 +337,14 @@ def check_c_res(): print("check B result") if mode_g: print(bb1 / (((loq + ab * nloq) @ Kqg + (ab * nlog) @ Kgg) @ f)) - plt.plot(bb1, label="b-fk") - plt.plot(((loq + ab * nloq) @ (Kqg) + (ab * nlog) @ (Kgg)) @ f, label="ana_pred") + # plt.plot(bb1, label="b-fk") + # plt.plot(((loq + ab * nloq) @ (Kqg) + (ab * nlog) @ (Kgg)) @ f, label="ana_pred") else: print(bb1 / (((loc + ab * nloc) @ Kqq + (ab * nlog) @ Kgq) @ f)) - plt.plot(bb1, label="b-fk") - plt.plot(((loc + ab * nloc) @ Kqq + (ab * nlog) @ Kgq) @ f, label="ana_pred") -plt.legend() -plt.show() + # plt.plot(bb1, label="b-fk") + # plt.plot(((loc + ab * nloc) @ Kqq + (ab * nlog) @ Kgq) @ f, label="ana_pred") +# plt.legend() +# plt.show() # check difference between B and C print("check B - C") if order_mode == Order.nlo or order_mode == Order.nnlo: @@ -330,26 +356,107 @@ def check_c_res(): if mode_g: # plt.plot(ab**2 * L * ((nloc @ pqg0 @ f)), label="pqg0") # plt.plot(ab**2 * L * ((nlog @ pgg0 @ f)), label="pgg0") - diff_ana = (ab * ab * L * ((nloq @ pqg0) + (nlog @ pgg0))) @ f - if signed: - diff_ana = ( - -2 * ab * L * (loq @ pqg0) - - ab * ab * L * ((nloq @ pqg0) + (nlog @ pgg0)) + if order_mode == Order.nlo: + diff_ana = (ab * ab * L * ((nloq @ pqg0) + (nlog @ pgg0))) @ f + if signed: + diff_ana = ( + -2 * ab * L * (loq @ pqg0) + - ab * ab * L * ((nloq @ pqg0) + (nlog @ pgg0)) + ) @ f + else: + nnlomix = ( + nlog @ pgg0 @ (pgg0 - beta0I) + + nlog @ pgq0 @ (pqg0 - beta0I) + + nloq @ pqq0 @ (pqg0 - beta0I) + + nloq @ pqg0 @ (pgg0 - beta0I) + ) + nnlomix2 = ( + nnlog @ pgg0 @ (pgg0 - beta0I) + + nnlog @ pgq0 @ (pqg0 - beta0I) + + nnloq @ pqq0 @ (pqg0 - beta0I) + + nnloq @ pqg0 @ (pgg0 - beta0I) + ) + diff_ana_nnlo = ( + ab + * ab + * ab + * ( + L + * ((nloq @ pqg1) + (nlog @ pgg1) + (nnlog @ pgg0) + (nnloq @ pqg0)) + + L * L * ((-1.0 / 2.0) * nnlomix) + ) + + ab + * ab + * ab + * ab + * ( + L * ((nnlog @ pgg1) + (nnloq @ pqg1)) + + L * L * ((-1.0 / 2.0) * (nnlomix2)) + ) ) @ f + diff_ana = diff_ana_nnlo + if signed: + diff_ana_nlo = ( + -2 * ab * L * (loq @ pqg0) + - ab * ab * L * ((nloq @ pqg0) + (nlog @ pgg0)) + ) @ f + diff_ana_nnlo = 0.0 + diff_ana = diff_ana_nlo + diff_ana_nnlo else: # plt.plot(a**2 * L * ((nloc @ pqq0 @ f)), label="pqq0") # plt.plot(a**2 * L * ((nlog @ pgq0 @ f)), label="pgq0") - diff_ana = (ab * ab * L * ((nloc @ pqq0) + (nlog @ pgq0))) @ f - if signed: - diff_ana = ( - -2 * ab * L * (loc @ pqq0) - - ab * ab * L * ((nloc @ pqq0) + (nlog @ pgq0)) + if order_mode == Order.nlo: + diff_ana = (ab * ab * L * ((nloc @ pqq0) + (nlog @ pgq0))) @ f + if signed: + diff_ana = ( + -2 * ab * L * (loc @ pqq0) + - ab * ab * L * ((nloc @ pqq0) + (nlog @ pgq0)) + ) @ f + else: + nnlomix = ( + nlog @ pgg0 @ (pgq0 - beta0I) + + nlog @ pgq0 @ (pqq0 - beta0I) + + nloq @ pqq0 @ (pqq0 - beta0I) + + nloq @ pqg0 @ (pgq0 - beta0I) + ) + nnlomix2 = ( + nnlog @ pgg0 @ (pgq0 - beta0I) + + nnlog @ pgq0 @ (pqq0 - beta0I) + + nnloq @ pqq0 @ (pqq0 - beta0I) + + nnloq @ pqg0 @ (pgq0 - beta0I) + ) + diff_ana_nnlo = ( + ab + * ab + * ab + * ( + L + * ((nloq @ pqq1) + (nlog @ pgq1) + (nnlog @ pgq0) + (nnloq @ pqq0)) + + L * L * ((-1.0 / 2.0) * nnlomix) + ) + + ab + * ab + * ab + * ab + * ( + L * ((nnlog @ pgq1) + (nnloq @ pqq1)) + + L * L * ((-1.0 / 2.0) * (nnlomix2)) + ) ) @ f + diff_ana = diff_ana_nnlo + if signed: + diff_ana_nlo = ( + -2 * ab * L * (loc @ pqq0) + - ab * ab * L * ((nloc @ pqq0) + (nlog @ pgq0)) + ) @ f + diff_ana_nnlo = 0.0 + diff_ana = diff_ana_nlo + diff_ana_nnlo + print((diff_grid / bb1) * 100) print(((diff_grid / diff_ana) - 1.0) * 100) plt.plot(diff_grid, label="grid") plt.plot(diff_ana, label="ana") plt.legend() plt.show() - plt.plot(((diff_grid[1:-5] / diff_ana[1:-5]) - 1.0) * 100, label="ratio") + plt.plot(((diff_grid / diff_ana) - 1.0) * 100, label="ratio") plt.legend() plt.show() From 36f147dda37b291a2d237bd2dcd58c727b5e67f3 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Thu, 10 Nov 2022 12:37:12 +0100 Subject: [PATCH 22/24] Add some comments --- debug3.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/debug3.py b/debug3.py index 31beed2d..3dff1f02 100644 --- a/debug3.py +++ b/debug3.py @@ -176,10 +176,10 @@ class Order(Enum): pgg0 = convolute_operator(nlo.pgg0(5), interp)[0].T pqq1 = convolute_operator(nlo.pqq1(5), interp)[0].T # pgg1 = convolute_operator(nlo.pgg1(5), interp)[0].T -pgg1 = pgg0 +pgg1 = pgg0 # just to have something, to be substituted with the correct one pqg1 = convolute_operator(nlo.pqg1(5), interp)[0].T # pgq1 = convolute_operator(nlo.pgq1(5), interp)[0].T -pgq1 = pgq0 +pgq1 = pgq0 # just to have something, to be substituted with the correct one beta0I = beta_qcd_as2(5) * np.eye(pqg0.shape[0]) beta0 = beta_qcd_as2(5) q2 = 300.0 From a55b2d773d58008527ba37a07d1cf084d0dac0e9 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Fri, 11 Nov 2022 14:19:49 +0100 Subject: [PATCH 23/24] Allow nnlo pdfs and different possible scheme B --- debug3.py | 48 ++++++++++++++++++++++++++++++++++-------------- 1 file changed, 34 insertions(+), 14 deletions(-) diff --git a/debug3.py b/debug3.py index 3dff1f02..d4bde103 100644 --- a/debug3.py +++ b/debug3.py @@ -52,10 +52,16 @@ class Order(Enum): signed = True # load respective PDF set -if mode_g: - pdf = lhapdf.mkPDF("gonly", 0) +if order_mode == Order.nnlo: + if mode_g: + pdf = lhapdf.mkPDF("gonly_nnlo", 0) + else: + pdf = lhapdf.mkPDF("conly_nnlo", 0) else: - pdf = lhapdf.mkPDF("conly", 0) + if mode_g: + pdf = lhapdf.mkPDF("gonly", 0) + else: + pdf = lhapdf.mkPDF("conly", 0) # load FK tables if order_mode == Order.nonlo: b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/no-nlo-test.pineappl.lz4") @@ -64,11 +70,23 @@ class Order(Enum): b1 = pineappl.fk_table.FkTable.read("data/fktables/2205/test.pineappl.lz4") c1 = pineappl.fk_table.FkTable.read("data/fktables/3205/test.pineappl.lz4") elif order_mode == Order.nnlo: + # 2405: bugged version + # 2505: the minus sign is globally in K + # 2605: the minus sign is in the definition of the gamma + # 2705: the minus sign is in the definition of L + # 2805: the minus sign is in the definition of beta0 and in the global K # b1 = pineappl.fk_table.FkTable.read("data/fktables/2405/test-nnlo.pineappl.lz4") #this is the one with the wrong sign c1 = pineappl.fk_table.FkTable.read("data/fktables/3405/test-nnlo.pineappl.lz4") - b1 = pineappl.fk_table.FkTable.read( - "data/fktables/2505/test-nnlo.pineappl.lz4" - ) # this is the one produced with the tentative solution of the sign + b1 = pineappl.fk_table.FkTable.read("data/fktables/2505/test-nnlo.pineappl.lz4") + # b1 = pineappl.fk_table.FkTable.read( + # "data/fktables/2605/test-nnlo.pineappl.lz4" + # ) + # b1 = pineappl.fk_table.FkTable.read( + # "data/fktables/2705/test-nnlo.pineappl.lz4" + # ) + # b1 = pineappl.fk_table.FkTable.read( + # "data/fktables/2805/test-nnlo.pineappl.lz4" + # ) out = yadism.output.Output.load_tar("../test0.tar") # compute predictions from FK bb1 = b1.convolute_with_one(2212, pdf.xfxQ2) @@ -181,7 +199,6 @@ class Order(Enum): # pgq1 = convolute_operator(nlo.pgq1(5), interp)[0].T pgq1 = pgq0 # just to have something, to be substituted with the correct one beta0I = beta_qcd_as2(5) * np.eye(pqg0.shape[0]) -beta0 = beta_qcd_as2(5) q2 = 300.0 muf2 = q2 * 2.0**2 # load theory card @@ -349,8 +366,8 @@ def check_c_res(): print("check B - C") if order_mode == Order.nlo or order_mode == Order.nnlo: diff_grid = bb1 - cc1 - plt.plot(bb1, label="b-grid") - plt.plot(cc1, label="c-grid") + plt.plot(bb1[3:], label="b-grid") + plt.plot(cc1[3:], label="c-grid") plt.legend() plt.show() if mode_g: @@ -364,6 +381,7 @@ def check_c_res(): - ab * ab * L * ((nloq @ pqg0) + (nlog @ pgg0)) ) @ f else: + pgg1 = pgg0 nnlomix = ( nlog @ pgg0 @ (pgg0 - beta0I) + nlog @ pgq0 @ (pqg0 - beta0I) @@ -413,6 +431,7 @@ def check_c_res(): - ab * ab * L * ((nloc @ pqq0) + (nlog @ pgq0)) ) @ f else: + pgq1 = pgq0 nnlomix = ( nlog @ pgg0 @ (pgq0 - beta0I) + nlog @ pgq0 @ (pqq0 - beta0I) @@ -451,12 +470,13 @@ def check_c_res(): ) @ f diff_ana_nnlo = 0.0 diff_ana = diff_ana_nlo + diff_ana_nnlo - print((diff_grid / bb1) * 100) - print(((diff_grid / diff_ana) - 1.0) * 100) - plt.plot(diff_grid, label="grid") - plt.plot(diff_ana, label="ana") + print((diff_grid[3:] / bb1[3:]) * 100) + print((diff_grid[3:] / cc1[3:]) * 100) + print(((diff_grid[3:] / diff_ana[3:]) - 1.0) * 100) + plt.plot(diff_grid[3:], label="grid") + plt.plot(diff_ana[3:], label="ana") plt.legend() plt.show() - plt.plot(((diff_grid / diff_ana) - 1.0) * 100, label="ratio") + plt.plot(((diff_grid[3:] / diff_ana[3:]) - 1.0) * 100, label="ratio") plt.legend() plt.show() From 99f6c78234c3a9d544ff0723038947a1e8bb2082 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Fri, 11 Nov 2022 14:23:06 +0100 Subject: [PATCH 24/24] Black --- debug3.py | 29 ++++++++++++----------------- 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/debug3.py b/debug3.py index d4bde103..86e71d6f 100644 --- a/debug3.py +++ b/debug3.py @@ -395,21 +395,21 @@ def check_c_res(): + nnloq @ pqg0 @ (pgg0 - beta0I) ) diff_ana_nnlo = ( - ab - * ab - * ab + (ab * ab * ab) * ( L * ((nloq @ pqg1) + (nlog @ pgg1) + (nnlog @ pgg0) + (nnloq @ pqg0)) + L * L * ((-1.0 / 2.0) * nnlomix) ) - + ab - * ab - * ab - * ab - * ( - L * ((nnlog @ pgg1) + (nnloq @ pqg1)) - + L * L * ((-1.0 / 2.0) * (nnlomix2)) + + ( + ab + * ab + * ab + * ab + * ( + L * ((nnlog @ pgg1) + (nnloq @ pqg1)) + + L * L * ((-1.0 / 2.0) * (nnlomix2)) + ) ) ) @ f diff_ana = diff_ana_nnlo @@ -445,18 +445,13 @@ def check_c_res(): + nnloq @ pqg0 @ (pgq0 - beta0I) ) diff_ana_nnlo = ( - ab - * ab - * ab + (ab * ab * ab) * ( L * ((nloq @ pqq1) + (nlog @ pgq1) + (nnlog @ pgq0) + (nnloq @ pqq0)) + L * L * ((-1.0 / 2.0) * nnlomix) ) - + ab - * ab - * ab - * ab + + (ab * ab * ab * ab) * ( L * ((nnlog @ pgq1) + (nnloq @ pqq1)) + L * L * ((-1.0 / 2.0) * (nnlomix2))