adding a hypercube transform function

enterprise/signals/gp_signals.py

@@ -13,9 +13,9 @@
 from sksparse.cholmod import cholesky
 
 from enterprise.signals import parameter, selections, signal_base, utils
+from enterprise.signals.signal_base import KernelMatrix
 from enterprise.signals.parameter import function
 from enterprise.signals.selections import Selection
-from enterprise.signals.utils import KernelMatrix
 
 # logging.basicConfig(format="%(levelname)s: %(name)s: %(message)s", level=logging.INFO)
 logger = logging.getLogger(__name__)
@@ -225,11 +225,11 @@ def get_timing_model_basis(use_svd=False, normed=True):
         return utils.unnormed_tm_basis()
 
 
-def TimingModel(coefficients=False, name="linear_timing_model", use_svd=False, normed=True):
+def TimingModel(coefficients=False, name="linear_timing_model", use_svd=False, normed=True, prior_variance=1e40):
     """Class factory for marginalized linear timing model signals."""
 
     basis = get_timing_model_basis(use_svd, normed)
-    prior = utils.tm_prior()
+    prior = utils.tm_prior(variance=prior_variance)
 
     BaseClass = BasisGP(prior, basis, coefficients=coefficients, name=name)
 
@@ -848,6 +848,7 @@ def MNMMNF(self, T):
 
     # we're ignoring logdet = True for two-dimensional cases, but OK
     def solve(self, right, left_array=None, logdet=False):
+        # compute generalized version of r+ N^-1 r
         if right.ndim == 1 and left_array is right:
             res = right
 
@@ -856,11 +857,13 @@ def solve(self, right, left_array=None, logdet=False):
             MNr = self.MNr(res)
             ret = rNr - np.dot(MNr, self.cf(MNr))
             return (ret, logdet_N + self.cf.logdet() + self.Mprior) if logdet else ret
+        # compute generalized version of T+ N^-1 r
         elif right.ndim == 1 and left_array is not None and left_array.ndim == 2:
             res, T = right, left_array
 
             TNr = self.Nmat.solve(res, left_array=T)
             return TNr - np.tensordot(self.MNMMNF(T), self.MNr(res), (0, 0))
+        # compute generalized version of T+ N^-1 T
         elif right.ndim == 2 and left_array is right:
             T = right
 

enterprise/signals/parameter.py

@@ -6,6 +6,7 @@
 
 import numpy as np
 from scipy.special import erf as _erf
+import scipy.stats as sstats
 
 from enterprise.signals.selections import selection_func
 
@@ -49,9 +50,15 @@ def __init__(self, name):
             msg = "Parameter classes need to define _prior, or _logprior."
             raise AttributeError(msg)
 
+        if hasattr(self, "_ppf"):
+            self.ppf = self._ppf(name)
+        else:
+            self.ppf = None
+
         self.type = self.__class__.__name__.lower()
 
     def get_logpdf(self, value=None, **kwargs):
+        # RvH: This exception cannot be triggered
         if not isinstance(self, Parameter):
             raise TypeError("You can only call get_logpdf() on an " "instantiated (named) Parameter.")
 
@@ -66,6 +73,7 @@ def get_logpdf(self, value=None, **kwargs):
         return logpdf if self._size is None else np.sum(logpdf)
 
     def get_pdf(self, value=None, **kwargs):
+        # RvH: This exception cannot be triggered
         if not isinstance(self, Parameter):
             raise TypeError("You can only call get_pdf() on an " "instantiated (named) Parameter.")
 
@@ -87,13 +95,23 @@ def sample(self, **kwargs):
             raise AttributeError("No sampler was provided for this Parameter.")
         else:
             if self.name in kwargs:
+                # RvH: This exception cannot be triggered
                 raise ValueError("You shouldn't give me my value when you're sampling me.!")
 
             if hasattr(self, "prior"):
                 return self.prior(func=self._sampler, size=self._size, **kwargs)
             else:
                 return self.logprior(func=self._sampler, size=self._size, **kwargs)
 
+    def get_ppf(self, value=None, **kwargs):
+        if self.ppf is None:
+            raise NotImplementedError("No ppf was implemented for this Parameter.")
+
+        if value is None and "params" in kwargs:
+            value = kwargs["params"][self.name]
+
+        return self.ppf(value, **kwargs)
+
     @property
     def size(self):
         return self._size
@@ -136,7 +154,7 @@ class GPCoefficients(Parameter):
     return GPCoefficients
 
 
-def UserParameter(prior=None, logprior=None, sampler=None, size=None):
+def UserParameter(prior=None, logprior=None, sampler=None, ppf=None, size=None):
     """Class factory for UserParameter, implementing Enterprise parameters
     with arbitrary priors. The prior is specified by way of an Enterprise
     ``Function`` of the form ``prior(value, [par1, par2])``. Optionally,
@@ -147,6 +165,7 @@ def UserParameter(prior=None, logprior=None, sampler=None, size=None):
     :param prior:   parameter prior pdf, given as Enterprise ``Function``
     :param sampler: function returning a randomly sampled parameter according
                     to prior
+    :param ppf:     percentage point function (inverse cdf), for this parameter
     :param size:    length for vector parameter
     :return:        ``UserParameter`` class
     """
@@ -157,6 +176,8 @@ class UserParameter(Parameter):
             _prior = prior
         if logprior is not None:
             _logprior = logprior
+        if ppf is not None:
+            _ppf = ppf
         _sampler = None if sampler is None else staticmethod(sampler)
         _typename = "UserParameter"
 
@@ -189,6 +210,12 @@ def UniformSampler(pmin, pmax, size=None):
     return np.random.uniform(pmin, pmax, size=size)
 
 
+def UniformPPF(value, pmin, pmax):
+    """Percentage Point function for Uniform paramters."""
+
+    return sstats.uniform.ppf(value, loc=pmin, scale=pmax - pmin)
+
+
 def Uniform(pmin, pmax, size=None):
     """Class factory for Uniform parameters (with pdf(x) ~ 1/[pmax - pmin]
     inside [pmin,pmax], 0 outside. Handles vectors correctly,
@@ -204,6 +231,7 @@ def Uniform(pmin, pmax, size=None):
     class Uniform(Parameter):
         _size = size
         _prior = Function(UniformPrior, pmin=pmin, pmax=pmax)
+        _ppf = Function(UniformPPF, pmin=pmin, pmax=pmax)
         _sampler = staticmethod(UniformSampler)
         _typename = _argrepr("Uniform", pmin=pmin, pmax=pmax)
 
@@ -233,6 +261,17 @@ def NormalSampler(mu, sigma, size=None):
         return np.random.normal(mu, sigma, size=size)
 
 
+def NormalPPF(value, mu, sigma):
+    """Prior function for Normal parameters.
+    Handles scalar mu and sigma, compatible vector value/mu/sigma,
+    vector value/mu and compatible covariance matrix sigma."""
+
+    if np.ndim(sigma) == 2:
+        raise NotImplementedError("PPF not implemented when sigma is 2D")
+
+    return sstats.norm.ppf(value, loc=mu, scale=sigma)
+
+
 def Normal(mu=0, sigma=1, size=None):
     """Class factory for Normal parameters (with pdf(x) ~ N(``mu``,``sigma``)).
     Handles vectors correctly if ``size == len(mu) == len(sigma)``,
@@ -249,6 +288,7 @@ def Normal(mu=0, sigma=1, size=None):
     class Normal(Parameter):
         _size = size
         _prior = Function(NormalPrior, mu=mu, sigma=sigma)
+        _ppf = Function(NormalPPF, mu=mu, sigma=sigma)
         _sampler = staticmethod(NormalSampler)
         _typename = _argrepr("Normal", mu=mu, sigma=sigma)
 
@@ -346,6 +386,13 @@ def LinearExpSampler(pmin, pmax, size=None):
     return np.log10(np.random.uniform(10**pmin, 10**pmax, size))
 
 
+def LinearExpPPF(value, pmin, pmax):
+    """Percentage Point function for Uniform paramters."""
+
+    ev = sstats.uniform.ppf(value, loc=10**pmin, scale=10**pmax - 10**pmin)
+    return np.log10(ev)
+
+
 def LinearExp(pmin, pmax, size=None):
     """Class factory for LinearExp parameters (with pdf(x) ~ 10^x,
     and 0 outside [``pmin``,``max``]). Handles vectors correctly
@@ -361,6 +408,7 @@ def LinearExp(pmin, pmax, size=None):
     class LinearExp(Parameter):
         _size = size
         _prior = Function(LinearExpPrior, pmin=pmin, pmax=pmax)
+        _ppf = Function(LinearExpPPF, pmin=pmin, pmax=pmax)
         _sampler = staticmethod(LinearExpSampler)
         _typename = _argrepr("LinearExp", pmin=pmin, pmax=pmax)
 
@@ -439,7 +487,6 @@ def __init__(self, name, psr=None):
 
             for kw, arg in self.func_kwargs.items():
                 if isinstance(arg, type) and issubclass(arg, (Parameter, ConstantParameter)):
-
                     # parameter name template:
                     #   pname_[signalname_][fname_]parname
                     pnames = [name, fname, kw]
@@ -582,7 +629,6 @@ def wrapper(*args, **kwargs):
                 and issubclass(arg, FunctionBase)
                 or isinstance(arg, FunctionBase)
             ):
-
                 return Function(func, **kwargs)
 
         # otherwise, we simply call the function

enterprise/signals/signal_base.py

@@ -24,7 +24,6 @@
 from enterprise.signals.parameter import Function  # noqa: F401
 from enterprise.signals.parameter import function  # noqa: F401
 from enterprise.signals.parameter import ConstantParameter
-from enterprise.signals.utils import KernelMatrix
 
 from enterprise import __version__
 from sys import version
@@ -677,6 +676,13 @@ def get_lnprior(self, params):
         params = params if isinstance(params, dict) else self.map_params(params)
 
         return np.sum([p.get_logpdf(params=params) for p in self.params])
+
+    def get_hypercube_transform(self, params):
+        # transform from unit cube to prior cube for nested sampling using PPFs
+        # map parameter vector if needed
+        params = params if isinstance(params, dict) else self.map_params(params)
+
+        return np.array([p.get_ppf(params=params) for p in self.params])
 
     @property
     def pulsars(self):
@@ -1212,6 +1218,89 @@ def solve(self, other, left_array=None, logdet=False):
         return (ret, self._get_logdet()) if logdet else ret
 
 
+class KernelMatrix(np.ndarray):
+    def __new__(cls, init):
+        if isinstance(init, int):
+            ret = np.zeros(init, "d").view(cls)
+        else:
+            ret = init.view(cls)
+
+        if ret.ndim == 2:
+            ret._cliques = -1 * np.ones(ret.shape[0])
+            ret._clcount = 0
+
+        return ret
+
+    # see PTA._setcliques
+    def _setcliques(self, idxs):
+        allidx = set(self._cliques[idxs])
+        maxidx = max(allidx)
+
+        if maxidx == -1:
+            self._cliques[idxs] = self._clcount
+            self._clcount = self._clcount + 1
+        else:
+            self._cliques[idxs] = maxidx
+            if len(allidx) > 1:
+                self._cliques[np.in1d(self._cliques, allidx)] = maxidx
+
+    def add(self, other, idx):
+        if other.ndim == 2 and self.ndim == 1:
+            self = KernelMatrix(np.diag(self))
+
+        if self.ndim == 1:
+            self[idx] += other
+        else:
+            if other.ndim == 1:
+                self[idx, idx] += other
+            else:
+                self._setcliques(idx)
+                idx = (idx, idx) if isinstance(idx, slice) else (idx[:, None], idx)
+                self[idx] += other
+
+        return self
+
+    def set(self, other, idx):
+        if other.ndim == 2 and self.ndim == 1:
+            self = KernelMatrix(np.diag(self))
+
+        if self.ndim == 1:
+            self[idx] = other
+        else:
+            if other.ndim == 1:
+                self[idx, idx] = other
+            else:
+                self._setcliques(idx)
+                idx = (idx, idx) if isinstance(idx, slice) else (idx[:, None], idx)
+                self[idx] = other
+
+        return self
+
+    def inv(self, logdet=False):
+        if self.ndim == 1:
+            inv = 1.0 / self
+
+            if logdet:
+                return inv, np.sum(np.log(self))
+            else:
+                return inv
+        else:
+            try:
+                cf = sl.cho_factor(self)
+                inv = sl.cho_solve(cf, np.identity(cf[0].shape[0]))
+                if logdet:
+                    ld = 2.0 * np.sum(np.log(np.diag(cf[0])))
+            except np.linalg.LinAlgError:
+                u, s, v = np.linalg.svd(self)
+                inv = np.dot(u / s, u.T)
+                if logdet:
+                    ld = np.sum(np.log(s))
+            if logdet:
+                return inv, ld
+            else:
+                return inv
+
+
 class ShermanMorrison(object):
     """Custom container class for Sherman-morrison array inversion."""