From fb5dc394724e531fd34102b67f34ed86123a60f5 Mon Sep 17 00:00:00 2001
From: achiefa <amedeochiefa@gmail.com>
Date: Thu, 30 May 2024 15:34:27 +0100
Subject: [PATCH] Implemented fit with spline
---
validphys2/src/validphys/config.py | 4 +
.../theorycovariance/construction.py | 92 ++++++++++++++++++-
2 files changed, 95 insertions(+), 1 deletion(-)
diff --git a/validphys2/src/validphys/config.py b/validphys2/src/validphys/config.py
index 9de191713b..954f9a1135 100644
--- a/validphys2/src/validphys/config.py
+++ b/validphys2/src/validphys/config.py
@@ -1126,6 +1126,10 @@ class CoreConfig(configparser.Config):
from validphys.theorycovariance.construction import thcov_HT_4
return thcov_HT_4
+ elif ht_version == 5:
+ from validphys.theorycovariance.construction import thcov_HT_5
+
+ return thcov_HT_5
else:
from validphys.theorycovariance.construction import covs_pt_prescrip
diff --git a/validphys2/src/validphys/theorycovariance/construction.py b/validphys2/src/validphys/theorycovariance/construction.py
index c2e6abebb7..58d1f7e601 100644
--- a/validphys2/src/validphys/theorycovariance/construction.py
+++ b/validphys2/src/validphys/theorycovariance/construction.py
@@ -372,7 +372,6 @@ def thcov_HT_4(combine_by_type_ht, ht_coeff_1, ht_coeff_2, ht_pt_prescription =
x = np.array([])
Q2 = np.array([])
y = np.array([])
- N = np.array([])
for exp in process_info.namelist[proc]:
# Locate position of the experiment
@@ -420,6 +419,97 @@ def thcov_HT_4(combine_by_type_ht, ht_coeff_1, ht_coeff_2, ht_pt_prescription =
return covmats
+def thcov_HT_5(combine_by_type_ht, H2_list, HL_list, reverse=False):
+ "Same as `thcov_HT` but implementing theory covariance method for each node of the spline."
+ process_info = combine_by_type_ht
+ running_index_tot = 0
+ start_proc_by_exp = defaultdict(list)
+ deltas = defaultdict(list)
+ included_proc = ["DIS NC"]
+ excluded_exp = {"DIS NC" : ["NMC_NC_NOTFIXED_DW_EM-F2"]}
+ included_exp = {}
+ for proc in included_proc:
+ aux = []
+ for exp in process_info.namelist[proc]:
+ if exp not in excluded_exp[proc]:
+ aux.append(exp)
+ included_exp[proc] = aux
+
+ # ABMP parametrisation
+ x_abmp = [0.0, 0.1, 0.3, 0.5, 0.7, 0.9, 1]
+
+ # Check that H2_list and HL_list have the same size as x
+ if (len(H2_list) != len(x_abmp)) or (len(HL_list) != len(x_abmp)):
+ raise ValueError(f"The size of HT parameters does not match the number of nodes in the spline.")
+
+ def wrapper_to_splines(i):
+ if not reverse:
+ shifted_H2_list = [0 for k in range(len(x_abmp))]
+ shifted_HL_list = [0 for k in range(len(x_abmp))]
+ shifted_H2_list[i] = H2_list[i]
+ shifted_HL_list[i] = HL_list[i]
+ else:
+ shifted_H2_list = H2_list.copy()
+ shifted_HL_list = HL_list.copy()
+ shifted_H2_list[i] = 0
+ shifted_HL_list[i] = 0
+
+ H_2 = scint.CubicSpline(x_abmp, shifted_H2_list)
+ H_L = scint.CubicSpline(x_abmp, shifted_HL_list)
+ H_2 = np.vectorize(H_2)
+ H_L = np.vectorize(H_L)
+ return H_2, H_L
+
+ for proc in process_info.namelist.keys():
+ running_index_proc = 0
+ x = np.array([])
+ Q2 = np.array([])
+ y = np.array([])
+
+ for exp in process_info.namelist[proc]:
+ # Locate position of the experiment
+ size = process_info.sizes[exp]
+ start_proc_by_exp[exp] = running_index_tot
+ running_index_tot += size
+ running_index_proc += size
+
+ # Compute shifts only for a subset of processes
+ if proc in included_proc and exp in included_exp[proc]:
+ #central = process_info.preds[proc][1][start_proc_by_exp[exp] : size] # Probably this is deprecated
+ x = process_info.data[proc].T[0][running_index_proc - size : running_index_proc]
+ Q2 = process_info.data[proc].T[1][running_index_proc - size : running_index_proc]
+ y = process_info.data[proc].T[2][running_index_proc - size : running_index_proc]
+
+ if "SIGMA" in exp:
+ N_2, N_L = compute_normalisation_by_experiment(exp, x, y, Q2)
+
+ elif "F2" in exp:
+ N_2 = np.ones(shape=x.shape)
+ N_L = np.zeros(shape=x.shape)
+
+ else:
+ raise ValueError(f"The normalisation for the observable is not known.")
+
+ # Loop over the parameter
+ for i in range(len(x_abmp)):
+ H_L, H_2 = wrapper_to_splines(i)
+ deltas[f"({i+1}+,0)"] += [N_2 * H_2(x) / Q2]
+ deltas[f"(0,{i+1}+)"] += [N_L * H_L(x) / Q2]
+
+
+ # Construct theory covmat
+ covmats = defaultdict(list)
+ for proc1 in included_proc:
+ for proc2 in included_proc:
+ for i, exp1 in enumerate(included_exp[proc1]):
+ for j, exp2 in enumerate(included_exp[proc2]):
+ s = np.zeros(shape=(deltas["(1+,0)"][i].size, deltas["(1+,0)"][j].size))
+ for par in deltas.keys():
+ s += np.outer(deltas[par][i], deltas[par][j])
+ start_locs = (start_proc_by_exp[exp1], start_proc_by_exp[exp2])
+ covmats[start_locs] = s
+ return covmats
+
def compute_normalisation_by_experiment(experiment_name, x, y, Q2):
N_2 = np.zeros(shape=y.shape)
N_L = np.zeros(shape=y.shape)
--
GitLab