[WIP] Branch for closure test chi2 table

validphys2/src/validphys/config.py

@@ -347,7 +347,13 @@ def produce_fitunderlyinglaw(self, fit):
         """
         with self.set_context(ns=self._curr_ns.new_child({"fit": fit})):
             _, datacuts = self.parse_from_("fit", "closuretest", write=False)
-        underlyinglaw = self.parse_pdf(datacuts["fakepdf"])
+        # underlyinglaw = self.parse_pdf(datacuts["fakepdf"])
+        fakepdf = datacuts["fakepdf"]
+
+        if isinstance(fakepdf, str):
+            underlyinglaw = self.parse_pdf(fakepdf)
+        else:
+            underlyinglaw = fakepdf
         return {"pdf": underlyinglaw}
 
     @element_of("hyperscans")

validphys2/src/validphys/core.py

@@ -42,6 +42,7 @@
 from validphys.hyperoptplot import HyperoptTrial
 from validphys.utils import experiments_to_dataset_inputs
 from validphys.lhapdfset import LHAPDFSet
+# from validphys.pineparser import pineappl_reader
 
 log = logging.getLogger(__name__)
 
@@ -513,7 +514,11 @@ def load(self):
 
         fktables = []
         for p in self.fkspecs:
-            fktable = p.load()
+            try:
+                fktable = p.load()
+            except Exception as e:
+                from validphys.pineparser import pineappl_reader
+                fktable = pineappl_reader(p)
             #IMPORTANT: We need to tell the python garbage collector to NOT free the
             #memory owned by the FKTable on garbage collection.
             #TODO: Do this automatically

validphys2/src/validphys/covmats.py

@@ -410,8 +410,23 @@ def sqrt_covmat(covariance_matrix):
                          f"{dimensions[1]}")
 
     sqrt_diags = np.sqrt(np.diag(covariance_matrix))
-    correlation_matrix = covariance_matrix / sqrt_diags[:, np.newaxis] / sqrt_diags
+    # There are some zero entries in sqrt_diags, this gives errors when finding the cholesky decomposition.
+
+    if np.any(sqrt_diags == 0):
+        # Use safe method if there are zeros
+        outer = np.outer(sqrt_diags, sqrt_diags)
+        correlation_matrix = np.divide(
+            covariance_matrix, outer, out=np.zeros_like(covariance_matrix), where=outer!=0
+        )
+    else:
+        # Normal division if all variances are nonzero
+        correlation_matrix = covariance_matrix / sqrt_diags[:, np.newaxis] / sqrt_diags
+
+    # Always fix the diagonal
+    np.fill_diagonal(correlation_matrix, 1.0)
+    # correlation_matrix = covariance_matrix / sqrt_diags[:, np.newaxis] / sqrt_diags
     decomp = la.cholesky(correlation_matrix)
+
     sqrt_matrix = (decomp * sqrt_diags).T
     return sqrt_matrix
 

validphys2/src/validphys/results.py

@@ -42,6 +42,7 @@
 
 from validphys.n3fit_data_utils import parse_simu_parameters_names_CF
 
+from validphys.pineparser import pineappl_reader
 
 log = logging.getLogger(__name__)
 
@@ -385,7 +386,13 @@ def dataset_bsm_factor(dataset, pdf, read_bsm_facs):
     if parsed_bsm_facs is None:
         # We want an array of ones that ndata x nrep
         # where ndata is the number of post cut datapoints
-        ndata = len(dataset.load().get_cv())
+        try:
+            ndata = len(dataset.load().get_cv())
+        except Exception:
+            fkspec = dataset.fkspecs[0]
+            fkdata = pineappl_reader(fkspec)
+            ndata = fkdata.ndata
+
         nrep = len(pdf)
         return np.ones((ndata, nrep))