Merge pull request #151 from pcubillos/multinest_new

Multinest new

mc3/plots/colors.py

@@ -133,7 +133,7 @@ def rainbow_text(ax, texts, fontsize, colors=None, loc='above'):
 
 class Theme():
     """A monochromatic color theme from given color"""
-    def __init__(self, color, alpha_light=0.15, alpha_dark=0.5):
+    def __init__(self, color, alpha_light=0.15, alpha_dark=0.7):
         """
         Parameters
         ----------

mc3/stats/stats.py

@@ -467,7 +467,7 @@ def cred_region(posterior=None, quantile=0.6827, pdf=None, xpdf=None):
     return pdf, xpdf, HPDmin
 
 
-class ppf_uniform(object):
+class ppf_uniform():
     """
     Percent-point function (PPF) for a uniform function between
     pmin and pmax.  Also known as inverse CDF or quantile function.
@@ -504,57 +504,75 @@ def __call__(self, u):
         return (self.pmax-self.pmin)*u + self.pmin
 
 
-class ppf_gaussian(object):
+class ppf_gaussian():
     """
-    Percent-point function (PPF) for a two-sided Gaussian function
+    Percent-point function (PPF) for a Gaussian distribution
+    (with potentially assymetric standard deviations)
     Also known as inverse CDF or quantile function.
 
     Parameters
     ----------
     loc: Float
         Center of the Gaussian function.
-    lo: Float
+    sigma_lo: Float
         Left-sided standard deviation (for values x < loc).
-    up: Float
+    sigma_up: Float
         Right-sided standard deviation (for values x > loc).
-
-    Returns
-    -------
-    ppf: Callable
-        The Gaussian's PPF.
+    pmin: Float
+        Left-sided domain boundary of the PPF.
+    pmax: Float
+        Right-sided domain boundary of the PPF.
 
     Examples
     --------
     >>> import mc3.stats as ms
-    >>> ppf_g = ms.ppf_gaussian(0.0, 1.0, 1.0)
+    >>> ppf_g = ms.ppf_gaussian(2.0, 1.0, 1.0)
     >>> # The domain of the output function is (0,1):
-    >>> print(ppf_g(1e-10), ppf_g(0.5), ppf_g(1.0-1e-10))
-    (-6.361340902404056, 0.0, 6.361340889697422)
+    >>> print(ppf_g(0.5))
+    2.0
+
+    >>> print(ppf_g(0.158))
+    0.9972883349734507
+
     >>> # Also works for np.array inputs:
-    >>> print(ppf_g(np.array([1e-10, 0.5, 1-1e-10])))
-    [-6.3613409   0.          6.36134089]
+    >>> print(ppf_g(np.array([0.025, 0.158, 0.5, 0.6])))
+    [0.04003602 0.99728833 2.         2.2533471 ]
     """
-    def __init__(self, loc, lo, up):
+    def __init__(self, loc, sigma_lo, sigma_up, pmin=-np.inf, pmax=np.inf):
         self.loc = loc
-        self.lo = lo
-        self.up = up
+        self.sigma_lo = sigma_lo
+        self.sigma_up = sigma_up
+        self.pmin = pmin
+        self.pmax = pmax
+        a = (self.pmin - self.loc) / self.sigma_lo
+        b = (self.pmax - self.loc) / self.sigma_up
+        self.rv_lo = ss.truncnorm(a, b, loc=loc, scale=sigma_lo)
+        if sigma_up != sigma_lo:
+            self.rv_up = ss.truncnorm(a, b, loc=loc, scale=sigma_up)
+        self.u_threshold = sigma_lo/(sigma_lo+sigma_up)
+        self._ufactor1 = 1.0 + sigma_up/sigma_lo
+        self._ufactor2 = 1.0 + sigma_lo/sigma_up
 
     def __call__(self, u):
-        if np.isscalar(u) and u < self.lo/(self.lo+self.up):
-            return ss.norm.ppf(0.5*u*(self.lo+self.up)/self.lo,
-                scale=self.lo, loc=self.loc)
-        elif np.isscalar(u):
-            return ss.norm.ppf(1-0.5*(1-u)*(1+self.lo/self.up),
-                scale=self.up, loc=self.loc)
-        # else:
+        if self.sigma_lo == self.sigma_up:
+            return self.rv_lo.ppf(u)
+
+        if np.isscalar(u):
+            if u < self.u_threshold:
+                return self.rv_lo.ppf(0.5*u*self._ufactor1)
+            return self.rv_up.ppf(1.0-0.5*(1-u)*self._ufactor2)
+
         icdf = np.empty_like(u)
-        left = u < self.lo/(self.lo+self.up)
-        icdf[ left] = ss.norm.ppf(0.5*u[left]*(1+self.up/self.lo),
-            scale=self.lo, loc=self.loc)
-        icdf[~left] = ss.norm.ppf(1-0.5*(1-u[~left])*(1+self.lo/self.up),
-            scale=self.up, loc=self.loc)
+        left = u < self.u_threshold
+        icdf[left] = self.rv_lo.ppf(0.5*u[left]*self._ufactor1)
+        icdf[~left] = self.rv_up.ppf(1.0-0.5*(1-u[~left])*self._ufactor2)
         return icdf
 
+    def draw(self, size):
+        u = np.random.uniform(size=size)
+        samples = self.__call__(u)
+        return samples
+
 
 def dwt_daub4(array, inverse=False):
     """
@@ -601,29 +619,36 @@ class Loglike(object):
     (sign of an invalid parameter set), return a large-negative
     log likelihood (to reject the sample).
     """
-    def __init__(self, data, uncert, func, params, indp, pstep):
+    def __init__(self, data, uncert, func, params, args, pstep):
         self.data = data
         self.uncert = uncert
         self.func = func
         self.params = params
-        self.indp = indp
+        self.args = args
         self.pstep = pstep
         self.ifree = pstep>0
         self.ishare = np.where(pstep<0)[0]
+
+        # Pre-calculate the part outside chi-square:
+        self._uncert_logl = -0.5*np.sum(np.log(2.0*np.pi*self.uncert**2.0))
+
+
     def __call__(self, params):
         self.params[self.ifree] = params
         for s in self.ishare:
             self.params[s] = self.params[-int(self.pstep[s])-1]
-        model = self.func(self.params, *self.indp)
+
+        model = self.func(self.params, *self.args)
         log_like = -0.5 * np.sum(
             ((self.data - model) / self.uncert)**2.0
         )
+        log_like += self._uncert_logl
         if not np.isfinite(log_like):
             log_like = -1.0e98
         return log_like
 
 
-class Prior_transform(object):
+class Prior_transform():
     """Wrapper to compute the PPF of a set of parameters."""
     def __init__(self, prior, priorlow, priorup, pmin, pmax, pstep):
         self.ppf = []
@@ -634,7 +659,7 @@ def __init__(self, prior, priorlow, priorup, pmin, pmax, pstep):
             if plo == 0.0 or pup == 0.0:
                 self.ppf.append(ppf_uniform(min, max))
             else:
-                self.ppf.append(ppf_gaussian(p0, plo, pup))
+                self.ppf.append(ppf_gaussian(p0, plo, pup, min, max))
     def __call__(self, u):
         return [ppf(v) for ppf,v in zip(self.ppf, u)]
 

mc3/version.py

@@ -4,6 +4,6 @@
 # mc3 Version:
 MC3_VER = 3  # Major version
 MC3_MIN = 1  # Minor version
-MC3_REV = 2  # Revision
+MC3_REV = 3  # Revision
 
 __version__ = f'{MC3_VER}.{MC3_MIN}.{MC3_REV}'

setup.py

@@ -22,7 +22,7 @@
 cfiles = list(filter(lambda x: not re.search('[.#].+[.]c$', x), cfiles))
 
 inc = [get_include(), incdir]
-eca = ['-ffast-math']
+eca = ['-O3', '-ffast-math']
 ela = []
 
 extensions = [

src_c/include/stats.h

@@ -36,12 +36,12 @@ Returns
 datarms: The root mean square of the data.
 ******************************************************************/
 double rms(double *data, const int n){
-  int i;
-  double datarms=0.0;
-  for (i=0; i<n; i++)
-    datarms += data[i]*data[i];
-  datarms = sqrt(datarms/n);
-  return datarms;
+    int i;
+    double datarms=0.0;
+    for (i=0; i<n; i++)
+        datarms += data[i]*data[i];
+    datarms = sqrt(datarms/n);
+    return datarms;
 }
 
 

tests/test_plots.py

@@ -110,10 +110,10 @@ def test_Theme():
     theme = mp.Theme(color)
 
     expected_light_color = np.array([0.25882353, 0.44705882, 0.90588235])
-    expected_dark_color = np.array([0.00588235, 0.13137255, 0.4372549 ])
+    expected_dark_color = np.array([0.00823529, 0.18392157, 0.61215686])
     expected_bad = expected_under = np.array([1., 1., 1., 1.])
     expected_first = np.array([0.85176471, 0.88941176, 0.98117647])
-    expected_last = np.array([0.00588235, 0.13137255, 0.4372549 ])
+    expected_last = np.array([0.00823529, 0.18392157, 0.61215686])
 
     assert theme.color == color
     np.testing.assert_allclose(theme.light_color, expected_light_color, rtol)