Merge pull request #1743 from NNPDF/hessian_pdf_covmat_for_theory_pre…

…dictions

Hessian PDF Covariance Matrix for theory predictions

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 @@
     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