Remove libNNPDF::LHAPDFSet for a vp-based class

doc/sphinx/source/vp/pydataobjs.rst

@@ -255,3 +255,36 @@ from the API::
         "theoryid": 162
     }
     total_cov = API.dataset_inputs_covmat_from_systematics(**inp)
+
+Loading LHAPDF PDFs
+-------------------
+
+A wrapper class for LHAPDF PDFs is implemented in the :py:mod:`validphys.lhapdfset` module.
+An instance of this module will provide with a handful of useful wrappers to the underlying
+LHAPDF python interface. This is also the output of the ``pdf.load()`` method.
+
+For example, the following will return the values for all 100 members of NNPDF4.0 for
+the gluon and the d-quark, at three values of ``x`` at ``Q=91.2``.
+
+.. code-block:: python
+
+    from validphys.api import API
+    pdf = API.pdf(pdf="NNPDF40_nnlo_as_01180")
+    l_pdf = pdf.load()
+    alpha_s = l_pdf.central_member.alphasQ(91.2)
+    results = l_pdf.grid_values([21,1], [0.1, 0.2, 0.3], [91.2])
+
+For Monte Carlo PDFs the result of the replica 0 (the average of all replicas) is not returned
+when using the ``grid_values`` method 
+(or in general when utilizing any action that computes a value per replica),
+for convenience, a ``enable_central_value`` is provided as part of the :py:class:`validphys.core.PDF`
+class:
+
+.. code-block:: python
+
+    from validphys.api import API
+    pdf = API.pdf(pdf="NNPDF40_nnlo_as_01180")
+    l_pdf = pdf.load()
+    with pdf.enable_central_value():
+        m101 = l_pdf.get_members()
+    m100 = l_pdf.get_members()

validphys2/src/validphys/app.py

@@ -14,7 +14,7 @@
 import logging
 import contextlib
 
-
+import lhapdf
 from reportengine import app
 
 from validphys.config import Config, Environment
@@ -129,7 +129,7 @@ def init(self):
         if not cout:
             import NNPDF
             NNPDF.SetVerbosity(0)
-            NNPDF.SetLHAPDFVerbosity(0)
+            lhapdf.setVerbosity(0)
 
     @staticmethod
     def upload_context(do_upload, output):

validphys2/src/validphys/core.py

@@ -9,6 +9,7 @@
 from __future__ import generator_stop
 
 from collections import namedtuple
+from contextlib import contextmanager
 import re
 import enum
 import functools
@@ -22,7 +23,7 @@
 from reportengine.baseexceptions import AsInputError
 from reportengine.compat import yaml
 
-from NNPDF import (LHAPDFSet,
+from NNPDF import (
     CommonData,
     FKTable,
     FKSet,
@@ -37,23 +38,10 @@
 from validphys.theorydbutils import fetch_theory
 from validphys.hyperoptplot import HyperoptTrial
 from validphys.utils import experiments_to_dataset_inputs
+from validphys.lhapdfset import LHAPDFSet
 
 log = logging.getLogger(__name__)
 
-
-#TODO: Remove this eventually
-#Bacward compatibility error type names
-#Swig renamed these for no reason whatsoever.
-try:
-    LHAPDFSet.erType_ER_EIG
-except AttributeError:
-    import warnings
-    warnings.warn("libnnpdf out of date. Setting backwards compatible names")
-    LHAPDFSet.erType_ER_MC = LHAPDFSet.ER_MC
-    LHAPDFSet.erType_ER_EIG = LHAPDFSet.ER_EIG
-    LHAPDFSet.erType_ER_EIG90 = LHAPDFSet.ER_EIG90
-    LHAPDFSet.erType_ER_SYMEIG = LHAPDFSet.ER_SYMEIG
-
 class TupleComp:
 
     @classmethod
@@ -91,10 +79,33 @@ def __dir__(self):
 PDFSETS = _PDFSETS()
 
 class PDF(TupleComp):
+    """Wrapper class for the validphys PDF object to easily manage
+    both Monte Carlo and Hessian sets
+
+    Offers a context manager (``enable_central_value``) to include the central value
+    in Monte Carlo sets.
+
+    Examples
+    --------
+    >>> from validphys.api import API
+    >>> from validphys.convolution import predictions
+    >>> args = {"dataset_input":{"dataset": "ATLASTTBARTOT"}, "theoryid":200, "use_cuts":"internal"}
+    >>> ds = API.dataset(**args)
+    >>> pdf = API.pdf(pdf="NNPDF40_nnlo_as_01180")
+    >>> with pdf.enable_central_value():
+    >>>     preds_with_cv = predictions(ds, pdf)
+    >>> preds_no_cv = predictions(ds, pdf)
+    >>> len(preds_with_cv.columns)
+    101
+    >>> len(preds_no_cv.columns)
+    100
+    """
 
     def __init__(self, name):
         self.name = name
         self._plotname = name
+        self._lhapdfset = None
+        self._include_cv = False
         super().__init__(name)
 
 
@@ -158,16 +169,16 @@ def rescale_factor(self):
         else:
             return 1
 
-    @functools.lru_cache(maxsize=16)
     def load(self):
-        return LHAPDFSet(self.name, self.nnpdf_error)
+        if self._lhapdfset is None:
+            self._lhapdfset = LHAPDFSet(self.name, self.nnpdf_error)
+        return self._lhapdfset
 
     @functools.lru_cache(maxsize=2)
     def load_t0(self):
-        """Load the PDF as a t0 set"""
+        """Reload the PDF as a t0 set"""
         return LHAPDFSet(self.name, LHAPDFSet.erType_ER_MCT0)
 
-
     def __str__(self):
         return self.label
 
@@ -178,7 +189,7 @@ def __len__(self):
 
     @property
     def nnpdf_error(self):
-        """Return the NNPDF error tag, used to build the `LHAPDFSet` objeect"""
+        """Return the NNPDF error tag, used to build the `LHAPDFSet` object"""
         error = self.ErrorType
         if error == "replicas":
             return LHAPDFSet.erType_ER_MC
@@ -217,6 +228,8 @@ def grid_values_index(self):
         len(pdf))`` for Monte Carlo sets, because replica 0 is not selected
         and ``range(0, len(pdf))`` for hessian sets.
 
+        If ``include_cv`` is set to True, add a central value column for member 0
+        for Monte Carlo error sets
 
         Returns
         -------
@@ -229,7 +242,10 @@ def grid_values_index(self):
         """
         err = self.nnpdf_error
         if err is LHAPDFSet.erType_ER_MC:
-            return range(1, len(self))
+            if self._include_cv:
+                return ["CV"] + list(range(1, len(self)))
+            else:
+                return range(1, len(self))
         elif err in (LHAPDFSet.erType_ER_SYMEIG, LHAPDFSet.erType_ER_EIG, LHAPDFSet.erType_ER_EIG90):
             return range(0, len(self))
         else:
@@ -242,6 +258,17 @@ def get_members(self):
         """
         return len(self.grid_values_index)
 
+    @contextmanager
+    def enable_central_value(self):
+        """Context manager within which the central value is included
+        regardless of the error type of the PDF set"""
+        # Get a reference to the base PDF set of this class
+        pdfset = self.load()
+        pdfset.include_cv = True
+        self._include_cv = True
+        yield
+        self._include_cv = False
+        pdfset.include_cv = False
 
 
 kinlabels_latex = CommonData.kinLabel_latex.asdict()

validphys2/src/validphys/covmats.py

@@ -489,7 +489,7 @@ def dataset_inputs_covmat(
         #Copy data to avoid chaos
         loaded_data = type(loaded_data)(loaded_data)
         log.debug("Setting T0 predictions for %s" % data)
-        loaded_data.SetT0(t0set.load_t0())
+        loaded_data.SetT0(t0set.load_t0().as_libNNPDF())
 
     covmat = loaded_data.get_covmat()
 
@@ -581,7 +581,7 @@ def covmat(
         #Copy data to avoid chaos
         loaded_data = type(loaded_data)(loaded_data)
         log.debug("Setting T0 predictions for %s" % dataset)
-        loaded_data.SetT0(t0set.load_t0())
+        loaded_data.SetT0(t0set.load_t0().as_libNNPDF())
 
     covmat = loaded_data.get_covmat()
     if fitthcovmat:

validphys2/src/validphys/filters.py

@@ -170,7 +170,7 @@ def _filter_closure_data(filter_path, data, fakepdfset, fakenoise, errorsize):
     # Load data, don't cache result
     loaded_data = data.load.__wrapped__(data)
     # generate level 1 shift if fakenoise
-    loaded_data.MakeClosure(fakeset, fakenoise)
+    loaded_data.MakeClosure(fakeset.as_libNNPDF(), fakenoise)
     for j, dataset in enumerate(data.datasets):
         path = filter_path / dataset.name
         nfull, ncut = _write_ds_cut_data(path, dataset)

validphys2/src/validphys/gridvalues.py

@@ -10,14 +10,8 @@
 
 import numpy as np
 
-from NNPDF import REALDOUBLE, LHAPDFSet
-
 from validphys.core import PDF
-
-#Numpy is unhappy about downcasting double to float implicitly, so we have
-#to manually make sure all input arrays correspond to NNPDF::real.
-FLTYPE = np.int32
-REALTYPE = np.double if REALDOUBLE else np.float32
+from validphys.lhapdfset import LHAPDFSet
 
 # Canonical ordering of PDG quark flavour codes
 QUARK_FLAVOURS = (-6, -5, -4, -3, -2, -1, 1, 2, 3, 4, 5, 6)
@@ -54,9 +48,9 @@
 
 def _grid_values(lpdf, flmat, xmat, qmat):
     """Compute lpdf.grid_values with more forgiving argument types"""
-    flmat = np.atleast_1d(np.asanyarray(flmat, dtype=FLTYPE))
-    xmat = np.atleast_1d(np.asarray(xmat, dtype=REALTYPE))
-    qmat = np.atleast_1d(np.asarray(qmat, dtype=REALTYPE))
+    flmat = np.atleast_1d(np.asanyarray(flmat))
+    xmat = np.atleast_1d(np.asarray(xmat))
+    qmat = np.atleast_1d(np.asarray(qmat))
     return lpdf.grid_values(flmat, xmat, qmat)
 
 def grid_values(pdf:PDF, flmat, xmat, qmat):