Apply comments

Author: achiefa

doc/sphinx/source/vp/theorycov/power_corrections.rst

@@ -155,7 +155,7 @@ The following keys are used:
   length of ``nodes``.
 - ``pc_included_procs``: list of process types to which power corrections
   are applied (e.g. ``["DIS NC", "DIS CC", "JETS", "DIJET"]``).
-- ``pc_excluded_exps``: list of dataset names to exclude from power corrections
+- ``pc_excluded_datasets``: list of dataset names to exclude from power corrections
   even if their process type is included.
 - ``pdf``: the PDF used for computing the theory predictions that enter the
   multiplicative shifts.
@@ -187,7 +187,7 @@ Example
           yshift: [2.0, 2.0, 2.0, 2.0, 2.0]
           nodes:  [0.25, 0.75, 1.25, 1.75, 2.25]
       pc_included_procs: ["JETS", "DIJET", "DIS NC", "DIS CC"]
-      pc_excluded_exps:
+      pc_excluded_datasets:
         - HERA_NC_318GEV_EAVG_CHARM-SIGMARED
         - HERA_NC_318GEV_EAVG_BOTTOM-SIGMARED
       pdf: NNPDF40_nnlo_as_01180
@@ -234,5 +234,5 @@ Module reference
 - ``covmat_power_corrections(deltas1, deltas2)``:
   computes the theory covariance sub-matrix between two datasets from their
   shift dictionaries.
-- ``covs_pt_prescrip_pc(combine_by_type, point_prescription, pdf, pc_parameters, pc_included_procs, pc_excluded_exps)``:
+- ``covs_pt_prescrip_pc(combine_by_type, point_prescription, pdf, pc_parameters, pc_included_procs, pc_excluded_datasets)``:
   assembles the full power correction covariance matrix across all datasets.

n3fit/runcards/examples/Basic_runcard_pc_covmat.yml

@@ -17,12 +17,11 @@ dataset_inputs:
 - {dataset: NUTEV_CC_NOTFIXED_FE_NU-SIGMARED, cfac: [MAS], variant: legacy_dw}
 - {dataset: NUTEV_CC_NOTFIXED_FE_NB-SIGMARED, cfac: [MAS], variant: legacy_dw}
 - {dataset: HERA_NC_318GEV_EM-SIGMARED}
-- {dataset: ATLAS_1JET_8TEV_R06_PTY, frac: 0.75, variant: decorrelated}
-- {dataset: ATLAS_2JET_7TEV_R06_M12Y, frac: 0.75}
-- {dataset: CMS_1JET_8TEV_PTY, frac: 0.75}
-- {dataset: CMS_2JET_7TEV_M12-Y, frac: 0.75}
-- {dataset: CMS_2JET_13TEV_M12-YSTAR-YB-R08, frac: 0.75}
-- {dataset: ATLAS_2JET_7TEV_R06_M12Y, frac: 0.75}
+- {dataset: ATLAS_1JET_8TEV_R06_PTY, variant: decorrelated}
+- {dataset: ATLAS_2JET_7TEV_R06_M12Y}
+- {dataset: CMS_1JET_8TEV_PTY}
+- {dataset: CMS_2JET_7TEV_M12-Y}
+- {dataset: CMS_2JET_13TEV_M12-YSTAR-YB-R08}
 
 ################################################################################
 diagonal_frac: 0.75
@@ -44,8 +43,8 @@ theorycovmatconfig:
     Hj:         {yshift: [2.0, 2.0, 2.0, 2.0, 2.0, 2.0], nodes: [0.25, 0.75, 1.25, 1.75, 2.25, 2.75]}
     H2j_ystar:  {yshift: [2.0, 2.0, 2.0, 2.0, 2.0, 2.0], nodes: [0.25, 0.75, 1.25, 1.75, 2.25, 2.75]}
     H2j_yb:   {yshift: [2.0, 2.0, 2.0, 2.0, 2.0], nodes: [0.25, 0.75, 1.25, 1.75, 2.25]}
-  pc_included_procs: ["DIJET"]
-  pc_excluded_exps: []
+  pc_included_procs: ["DIS NC", "DIS CC", "JETS", "DIJET"]
+  pc_excluded_datasets: []
   pdf: 260202-jk-nnpdf41-mhou
   use_thcovmat_in_fitting: true
   use_thcovmat_in_sampling: true
@@ -56,8 +55,8 @@ nnseed: 953262798
 mcseed: 1437981271
 genrep: true
 parameters: # This defines the parameter dictionary that is passed to the Model Trainer
-  nodes_per_layer: [70, 50, 25, 20, 9]
-  activation_per_layer: [tanh, tanh, tanh, tanh, linear]
+  nodes_per_layer: [25, 20, 9]
+  activation_per_layer: [tanh, tanh, linear]
   initializer: glorot_normal
   optimizer:
     clipnorm: 6.073e-6

validphys2/src/validphys/theorycovariance/construction.py

@@ -12,8 +12,6 @@
 
 from reportengine import collect
 from reportengine.table import table
-
-pass
 from validphys.checks import check_pc_parameters
 from validphys.results import results, results_central
 from validphys.theorycovariance.higher_twist_functions import compute_deltas_pc
@@ -78,7 +76,7 @@ def combine_by_type(each_dataset_results_central_bytheory, groups_data_by_proces
     dataset_size = defaultdict(list)
     theories_by_process = defaultdict(list)
     ordered_names = defaultdict(list)
-    data_spec = defaultdict(list)
+    data_spec = {}
     for dataset in each_dataset_results_central_bytheory:
         name = dataset[0][0].name
         theory_centrals = [x[1].central_value for x in dataset]
@@ -91,8 +89,8 @@ def combine_by_type(each_dataset_results_central_bytheory, groups_data_by_proces
 
     # Store DataGroupSpecs instances
     for group_proc in groups_data_by_process:
-        for exp_set in group_proc.datasets:
-            data_spec[exp_set.name] = exp_set
+        for ds in group_proc.datasets:
+            data_spec[ds.name] = ds
 
     process_info = ProcessInfo(
         preds=theories_by_process, namelist=ordered_names, sizes=dataset_size, data_spec=data_spec
@@ -362,7 +360,7 @@ def covs_pt_prescrip_mhou(combine_by_type, point_prescription):
 
 @check_pc_parameters
 def covs_pt_prescrip_pc(
-    combine_by_type, point_prescription, pdf, pc_parameters, pc_included_procs, pc_excluded_exps
+    combine_by_type, point_prescription, pdf, pc_parameters, pc_included_procs, pc_excluded_datasets
 ):
     """Produces the sub-matrices of the theory covariance matrix for power
     corrections. Sub-matrices correspond to applying power corrected shifts
@@ -372,33 +370,34 @@ def covs_pt_prescrip_pc(
     running_index = 0
 
     covmats = defaultdict(list)
-    start_proc_by_exp = defaultdict(list)
-    for exp_name, data_spec in datagroup_spec.items():
-        start_proc_by_exp[exp_name] = running_index
-        running_index += data_spec.load_commondata().ndata  # takes cuts into account
-
-    for exp_name1, data_spec1 in datagroup_spec.items():
-        for exp_name2, data_spec2 in datagroup_spec.items():
-            process_type1 = process_lookup(exp_name1)
-            process_type2 = process_lookup(exp_name2)
-
-            is_excluded_exp = any(name in pc_excluded_exps for name in [exp_name1, exp_name2])
-            is_included_proc = any(
+    start_proc_by_dsname = {}
+    for ds_name, ds_spec in datagroup_spec.items():
+        start_proc_by_dsname[ds_name] = running_index
+        running_index += ds_spec.load_commondata().ndata  # takes cuts into account
+
+    for ds_name1, ds_spec1 in datagroup_spec.items():
+        for ds_name2, ds_spec2 in datagroup_spec.items():
+            process_type1 = process_lookup(ds_name1)
+            process_type2 = process_lookup(ds_name2)
+
+            is_excluded_ds = any(name in pc_excluded_datasets for name in [ds_name1, ds_name2])
+            is_excluded_proc = any(
                 proc not in pc_included_procs for proc in [process_type1, process_type2]
             )
-            if not (is_excluded_exp or is_included_proc):
-                deltas1 = compute_deltas_pc(data_spec1, pdf, pc_parameters)
-                deltas2 = compute_deltas_pc(data_spec2, pdf, pc_parameters)
+            if is_excluded_ds or is_excluded_proc:
+                continue
 
-                # Convert deltas1 and deltas2 to the format expected by compute_covs_pt_prescrip
-                deltas1 = [np.array(deltas1[shift]) for shift in sorted(deltas1.keys())]
-                deltas2 = [np.array(deltas2[shift]) for shift in sorted(deltas2.keys())]
+            deltas1_dict = compute_deltas_pc(ds_spec1, pdf, pc_parameters)
+            deltas2_dict = compute_deltas_pc(ds_spec2, pdf, pc_parameters)
+
+            # Convert deltas1 and deltas2 to the format expected by compute_covs_pt_prescrip
+            deltas1 = [np.array(deltas1_dict[shift]) for shift in sorted(deltas1_dict.keys())]
+            deltas2 = [np.array(deltas2_dict[shift]) for shift in sorted(deltas2_dict.keys())]
+
+            s = compute_covs_pt_prescrip(point_prescription, ds_name1, deltas1, ds_name2, deltas2)
+            start_locs = (start_proc_by_dsname[ds_name1], start_proc_by_dsname[ds_name2])
+            covmats[start_locs] = s
 
-                s = compute_covs_pt_prescrip(
-                    point_prescription, exp_name1, deltas1, exp_name2, deltas2
-                )
-                start_locs = (start_proc_by_exp[exp_name1], start_proc_by_exp[exp_name2])
-                covmats[start_locs] = s
     return covmats
 
 

validphys2/src/validphys/theorycovariance/higher_twist_functions.py

@@ -118,8 +118,7 @@ def get_pc_type(
             )
         return ("H2j_ystar", "H2j_yb")
 
-    else:
-        raise RuntimeError(f"{process_type} has not been implemented.")
+    raise RuntimeError(f"{exp_name}: {process_type} has not been implemented.")
 
 
 def _get_dis_pc_type(exp_name: str) -> PCTypeResult:
@@ -139,23 +138,23 @@ def _get_dis_pc_type(exp_name: str) -> PCTypeResult:
     raise ValueError(f"The DIS observable for {exp_name} has not been implemented.")
 
 
-def _get_dijet_pc_type(experiment: str, pc_dict: Optional[dict] = None) -> str:
+def _get_dijet_pc_type(experiment: str, pc_dict: Optional[dict] = {}) -> str:
     """Resolve DIJET experiment to PC type key with optional fallback."""
     if pc_dict.get("H2j_ystar") and pc_dict.get("H2j_yb"):
         return ("H2j_ystar", "H2j_yb")
+
+    if experiment == 'ATLAS':
+        specific_key = "H2j_ystar"
+    elif experiment == 'CMS':
+        specific_key = "H2j_ymax"
     else:
-        if experiment == 'ATLAS':
-            specific_key = "H2j_ystar"
-        elif experiment == 'CMS':
-            specific_key = "H2j_ymax"
-        else:
-            raise ValueError(f"{experiment} is not implemented for DIJET.")
+        raise ValueError(f"{experiment} is not implemented for DIJET.")
 
     return specific_key
 
 
 def linear_bin_function(
-    a: npt.ArrayLike, y_shift: npt.ArrayLike, bin_edges: npt.ArrayLike
+    a: npt.ArrayLike, y_shift: npt.ArrayLike, nodes: npt.ArrayLike
 ) -> np.ndarray:
     """
     This function defines the linear bin function used to construct the prior. Specifically,
@@ -168,7 +167,7 @@ def linear_bin_function(
       A one-dimensional array of points at which the function is evaluated.
     y_shift:  ArrayLike of float
       A one-dimensional array whose elements represent the y-value of each bin
-    bin_nodes: ArrayLike of float
+    nodes: ArrayLike of float
       A one-dimensional array containing the edges of the bins. The bins are
       constructed using pairs of consecutive points.
 
@@ -177,36 +176,9 @@ def linear_bin_function(
     A one-dimensional array containing the function values evaluated at the points
     specified in `a`.
     """
-    res = np.zeros_like(a)
-    for shift_pos, shift in enumerate(y_shift):
-        if shift_pos > 0 and shift_pos < len(y_shift) - 1:
-            bin_low = bin_edges[shift_pos - 1]
-            bin_high = bin_edges[shift_pos + 1]
-            bin_mid = bin_edges[shift_pos]
-            m1 = shift / (bin_mid - bin_low)
-            m2 = shift / (bin_high - bin_mid)
-        elif shift_pos == 0:  # Left-most bin
-            bin_high = bin_edges[shift_pos + 1]
-            bin_mid = bin_edges[shift_pos]
-            bin_low = bin_mid
-            m1 = 0.0
-            m2 = shift / (bin_high - bin_mid)
-        else:  # Right-most bin
-            bin_low = bin_edges[shift_pos - 1]
-            bin_mid = bin_edges[shift_pos]
-            bin_high = bin_mid
-            m1 = shift / (bin_mid - bin_low)
-            m2 = 0.0
-        cond_low = np.multiply(
-            a >= bin_low, a < bin_mid if shift_pos != len(y_shift) - 1 else a <= bin_mid
-        )
-        cond_high = np.multiply(
-            a >= bin_mid, a < bin_high if shift_pos != len(y_shift) - 1 else a <= bin_high
-        )
-        res = np.add(res, [m1 * (val - bin_low) if cond else 0.0 for val, cond in zip(a, cond_low)])
-        res = np.add(
-            res, [-m2 * (val - bin_high) if cond else 0.0 for val, cond in zip(a, cond_high)]
-        )
+    a = np.asarray(a, dtype=float)
+    res = np.interp(a, nodes, y_shift)
+    res[(a < nodes[0]) | (a > nodes[-1])] = 0.0
     return res
 
 
@@ -256,10 +228,10 @@ def jets_pc_func(
       One-dimensional array containing the shifts for each bin.
     nodes: ArrayLike
      One-dimensional array containing the edges of the bins in rapidity.
-    rap: ArrayLike
-      List of rapidity points at which the power correction is evaluated.
     pT: ArrayLike
       List of pT points at which the power correction is evaluated.
+    rap: ArrayLike
+      List of rapidity points at which the power correction is evaluated.
 
     Returns
     -------
@@ -442,6 +414,8 @@ def _F(tx, ty):
 
 # ---------------------------------------------------------------------------
 # Power correction callable factories
+# If these functions are used to compute c-factors,
+# no need to multiply by the observable.
 # ---------------------------------------------------------------------------
 def mult_dis_pc(nodes, x, q2, dataset_sp, pdf, cfactors: bool = False) -> callable:
     """
@@ -526,8 +500,6 @@ def func(y_values_p, y_values_d):
         h2d = dis_pc_func(y_values_d, d_nodes, x, q2)
         h2p = dis_pc_func(y_values_p, p_nodes, x, q2)
         if cfactors:
-            # If this function is used to compute c-factors,
-            # no need to multiply by the observable
             num = np.sum([np.ones_like(h2d), h2d], axis=0)
             denom = np.sum([np.ones_like(h2p), h2p], axis=0)
             return (num, denom)
@@ -583,10 +555,20 @@ def mult_dijet_pc(
     fk = fkspec.load_with_cuts(cuts)
     th_preds = central_fk_predictions(fk, pdf)
 
+    if distrb == '3D':
+        if len(ystar) == 0 or len(yb) == 0:
+            raise ValueError("For 3D distributions, ystar and yb must be provided.")
+    elif distrb == "ystar":
+        if len(ystar) == 0:
+            raise ValueError("For ystar distributions, ystar must be provided.")
+    elif distrb == "ymax":
+        if len(ymax) == 0:
+            raise ValueError("For ymax distributions, ymax must be provided.")
+    else:
+        raise ValueError(f"Unknown distrb value: {distrb!r}. Expected '3D', 'ystar', or 'ymax'.")
+
     def func(ystar_shifts, y_b_shifts):
         if distrb == '3D':
-            if len(ystar) == 0 or len(yb) == 0:
-                raise ValueError("For 3D distributions, ystar and yb must be provided.")
             result = dijet3D_pc_func(
                 ystar_shifts=ystar_shifts,
                 yb_shifts=y_b_shifts,
@@ -597,8 +579,6 @@ def func(ystar_shifts, y_b_shifts):
                 m_jj=m_jj,
             )
         elif distrb == "ystar":
-            if len(ystar) == 0:
-                raise ValueError("For ystar distributions, ystar must be provided.")
             result = dijet_ystar_pc_func(
                 ystar_shifts=ystar_shifts,
                 yb_shifts=y_b_shifts,
@@ -608,8 +588,6 @@ def func(ystar_shifts, y_b_shifts):
                 ystar=ystar,
             )
         elif distrb == "ymax":
-            if len(ymax) == 0:
-                raise ValueError("For ymax distributions, ymax must be provided.")
             result = dijet_ymax_pc_func(
                 ystar_shifts=ystar_shifts,
                 yb_shifts=y_b_shifts,
@@ -618,10 +596,6 @@ def func(ystar_shifts, y_b_shifts):
                 m_jj=m_jj,
                 ymax=ymax,
             )
-        else:
-            raise ValueError(
-                f"Unknown distrb value: {distrb!r}. Expected '3D', 'ystar', or 'ymax'."
-            )
 
         if cfactors:
             return np.sum([np.ones_like(result), result], axis=0)
@@ -756,9 +730,11 @@ def compute_deltas_pc(
     pc_type = get_pc_type(exp_name, process_type, experiment=experiment, pc_dict=pc_dict)
 
     if process_type.startswith('DIS'):
-
-        x = dataset_sp.commondata.metadata.load_kinematics()['x'].to_numpy().reshape(-1)[cuts]
-        q2 = dataset_sp.commondata.metadata.load_kinematics()['Q2'].to_numpy().reshape(-1)[cuts]
+        try:
+            x = dataset_sp.commondata.metadata.load_kinematics()['x'].to_numpy().reshape(-1)[cuts]
+            q2 = dataset_sp.commondata.metadata.load_kinematics()['Q2'].to_numpy().reshape(-1)[cuts]
+        except KeyError as e:
+            raise ValueError(f"Required kinematic variable not found: {e}")
 
         # NMC ratio: special case with two PC types (proton and deuteron)
         if isinstance(pc_type, tuple):
@@ -791,7 +767,6 @@ def compute_deltas_pc(
             raise ValueError(f"{experiment} is not implemented for DIJET.")
 
         kinematics = dataset_sp.commondata.metadata.load_kinematics()
-        # TODO Excpetion rule for CMS_2JET_7TEV_M12-Y should be removed
         rap = kinematics[rap_var].to_numpy().reshape(-1)[cuts]
         m_jj = kinematics['m_jj'].to_numpy().reshape(-1)[cuts]
 
@@ -814,6 +789,8 @@ def compute_deltas_pc(
             # Both ystar and yb parameters contribute for all 2D distributions
             deltas.update(_apply_pars_combs(pars_combs, ["H2j_ystar", "H2j_yb"], pc_func))
 
+        # Fallback to the case where power corrections have same structure
+        # as single-inclusice jets, with pT replaced by m_jj.
         else:
             nodes = pc_dict[pc_type]['nodes']
             pc_func = mult_jet_pc(nodes, m_jj, rap, dataset_sp, pdf, cfactors=cfactors)