Merge branch 'master' into fix_evolven3fit_for_nf=3

n3fit/src/evolven3fit_new/eko_utils.py

@@ -22,7 +22,7 @@ EVOLVEN3FIT_CONFIGS_DEFAULTS_TRN = {
 }
 
 EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA = {
-    "ev_op_iterations": 10,
+    "ev_op_iterations": 30,
     "ev_op_max_order": (1, 0),
     "evolution_method": "iterate-exact",
     "inversion_method": "exact",

validphys2/src/validphys/checks.py

@@ -43,6 +43,16 @@ def check_pdf_is_montecarlo(ns, **kwargs):
         raise CheckError(f"Error type of PDF {pdf} must be 'replicas' and not {etype}")
 
 
+@make_argcheck
+def check_pdf_is_montecarlo_or_symmhessian(ns, **kwargs):
+    pdf = ns['pdf']
+    etype = pdf.error_type
+    check(
+        etype in {'replicas', 'symhessian'},
+        f"Error type of PDF {pdf} must be either 'replicas' or 'symmhessian' and not {etype}",
+    )
+
+
 @make_check
 def check_know_errors(ns, **kwargs):
     pdf = ns['pdf']

validphys2/src/validphys/covmats.py

@@ -15,7 +15,7 @@ from validphys.checks import (
     check_data_cuts_match_theorycovmat,
     check_dataset_cuts_match_theorycovmat,
     check_norm_threshold,
-    check_pdf_is_montecarlo,
+    check_pdf_is_montecarlo_or_symmhessian,
     check_speclabels_different,
 )
 from validphys.commondata import loaded_commondata_with_cuts
@@ -677,11 +677,15 @@ def groups_corrmat(groups_covmat):
     return mat
 
 
-@check_pdf_is_montecarlo
+@check_pdf_is_montecarlo_or_symmhessian
 def pdferr_plus_covmat(dataset, pdf, covmat_t0_considered):
     """For a given `dataset`, returns the sum of the covariance matrix given by
-    `covmat_t0_considered` and the PDF error: a covariance matrix estimated from the
-    replica theory predictions from a given monte carlo `pdf`
+    `covmat_t0_considered` and the PDF error:
+    - If the PDF error_type is 'replicas', a covariance matrix is estimated from
+      the replica theory predictions
+    - If the PDF error_type is 'symmhessian', a covariance matrix is estimated using
+      formulas from (mc2hessian) https://arxiv.org/pdf/1505.06736.pdf
+
 
     Parameters
     ----------
@@ -716,7 +720,22 @@ def pdferr_plus_covmat(dataset, pdf, covmat_t0_considered):
     True
     """
     th = ThPredictionsResult.from_convolution(pdf, dataset)
-    pdf_cov = np.cov(th.error_members, rowvar=True)
+
+    if pdf.error_type == 'replicas':
+        pdf_cov = np.cov(th.error_members, rowvar=True)
+
+    elif pdf.error_type == 'symmhessian':
+        rescale_fac = pdf._rescale_factor()
+        hessian_eigenvectors = th.error_members
+        central_predictions = th.central_value
+
+        # need to subtract the central set which is not the same as the average of the
+        # Hessian eigenvectors.
+        X = hessian_eigenvectors - central_predictions.reshape((central_predictions.shape[0], 1))
+        # need to rescale the Hessian eigenvectors in case the eigenvector confidence interval is not 68%
+        X = X / rescale_fac
+        pdf_cov = X @ X.T
+
     return pdf_cov + covmat_t0_considered
 
 

validphys2/src/validphys/photon/compute.py

@@ -17,7 +17,6 @@ from .constants import ED2, EU2, NC, NL
 
 log = logging.getLogger(__name__)
 
-Q_IN = 100
 FIATLUX_DEFAULT = {
     "apfel": False,
     "eps_base": 1e-5,  # precision on final integration of double integral.
@@ -71,6 +70,9 @@ class Photon:
         path_to_F2 = theoryid.path / "fastkernel/fiatlux_dis_F2.pineappl.lz4"
         path_to_FL = theoryid.path / "fastkernel/fiatlux_dis_FL.pineappl.lz4"
         self.path_to_eko_photon = theoryid.path / "eko_photon.tar"
+        with EKO.read(self.path_to_eko_photon) as eko:
+            self.q_in = np.sqrt(eko.mu20)
+        
 
         # set fiatlux
         self.lux = {}
@@ -129,7 +131,7 @@ class Photon:
         """
         # Compute photon PDF
         log.info(f"Computing photon")
-        photon_qin = np.array([self.lux[replica].EvaluatePhoton(x, Q_IN**2).total for x in XGRID])
+        photon_qin = np.array([self.lux[replica].EvaluatePhoton(x, self.q_in**2).total for x in XGRID])
         photon_qin += self.generate_errors(replica)
         # fiatlux computes x * gamma(x)
         photon_qin /= XGRID
@@ -141,8 +143,8 @@ class Photon:
             # it has to be done inside vp-setupfit
 
             # construct PDFs
-            pdfs_init = np.zeros((len(eko.rotations.inputpids), len(XGRID)))
-            for j, pid in enumerate(eko.rotations.inputpids):
+            pdfs_init = np.zeros((len(eko.bases.inputpids), len(XGRID)))
+            for j, pid in enumerate(eko.bases.inputpids):
                 if pid == 22:
                     pdfs_init[j] = photon_qin
                     ph_id = j
@@ -150,12 +152,11 @@ class Photon:
                     if pid not in self.luxpdfset.flavors:
                         continue
                     pdfs_init[j] = np.array(
-                        [self.luxpdfset.xfxQ(x, Q_IN, replica, pid) / x for x in XGRID]
+                        [self.luxpdfset.xfxQ(x, self.q_in, replica, pid) / x for x in XGRID]
                     )
 
             # Apply EKO to PDFs
-            q2 = eko.mu2grid[0]
-            with eko.operator(q2) as elem:
+            for _, elem in eko.items():
                 pdfs_final = np.einsum("ajbk,bk", elem.operator, pdfs_init)
 
         photon_Q0 = pdfs_final[ph_id]
@@ -188,7 +189,7 @@ class Photon:
         if not self.additional_errors:
             return None
         extra_set = self.additional_errors.load()
-        qs = [Q_IN] * len(XGRID)
+        qs = [self.q_in] * len(XGRID)
         res_central = np.array(extra_set.central_member.xfxQ(22, XGRID, qs))
         res = []
         for idx_member in range(101, 107 + 1):

validphys2/src/validphys/tests/photon/test_compute.py

@@ -8,6 +8,7 @@ from validphys.photon.compute import Photon, Alpha
 from validphys.core import PDF as PDFset
 
 from ..conftest import PDF
+from eko.io import EKO
 
 
 class FakeTheory:
@@ -68,6 +69,17 @@ class FakeFiatlux:
     def EvaluatePhoton(self, x, q):
         return self.res
 
+class FakeEKO:
+    def __init__(self, path):
+        self.path = path
+        self.mu20 = 100**2
+    
+    def __enter__(self):
+        return self
+    
+    def __exit__(self, exc_type: type, _exc_value, _traceback):
+        pass
+
 
 class FakeStructureFunction:
     def __init__(self, path, pdfs):
@@ -95,6 +107,7 @@ def test_parameters_init(monkeypatch):
     monkeypatch.setattr(structure_functions, "F2LO", FakeF2LO)
     monkeypatch.setattr(fiatlux, "FiatLux", FakeFiatlux)
     monkeypatch.setattr(Photon, "compute_photon_array", lambda *args: np.zeros(196))
+    monkeypatch.setattr(EKO, "read", FakeEKO)
 
     photon = Photon(FakeTheory(), fiatlux_runcard, [1, 2, 3])
     alpha = Alpha(FakeTheory().get_description())

validphys2/src/validphys/tests/test_regressions.py

@@ -67,22 +67,32 @@ def test_mcreplica(data_config):
 def test_expcovmat(data_config):
     return API.groups_covmat_no_table(**data_config)
 
+
 @make_table_comp(parse_exp_mat)
 def test_t0covmat(data_witht0_internal_cuts_config):
     return API.groups_covmat_no_table(**data_witht0_internal_cuts_config)
 
+
 @make_table_comp(parse_exp_mat)
 def test_expsqrtcovmat(data_config):
     return API.groups_sqrtcovmat(**data_config)
 
+
 @make_table_comp(parse_exp_mat)
 def test_t0sqrtcovmat(data_witht0_internal_cuts_config):
     return API.groups_sqrtcovmat(**data_witht0_internal_cuts_config)
 
+
 @make_table_comp(parse_exp_mat)
 def test_pdf_plus_exp_covmat(data_internal_cuts_config):
     return API.groups_covmat_no_table(use_pdferr=True, **data_internal_cuts_config)
 
+
+@make_table_comp(parse_exp_mat)
+def test_hessian_pdf_plus_exp_covmat(hessian_data_internal_cuts_config):
+    return API.groups_covmat_no_table(use_pdferr=True, **hessian_data_internal_cuts_config)
+
+
 @make_table_comp(sane_load)
 def test_predictions(data_internal_cuts_config):
     # TODO: ideally we would change the baseline to just be corresponding columns