v0.9.2

This commit fixes a number of bugs and also introduces a significant
amount of stability fixes for ellipsoidal decomposition to improve behavior
in the sparsely-sampled regime and/or the high-dimensional regime. All of
these changes have been stress-tested using 50-D to 250-D Gaussians, which have
now been added to the examples. A summary of changes is below:

Ellipsoids:
- Adds in support for computing shrunk covariances from `sklearn.covariance`
  using the OAS estimator (which is optimal for Gaussian distributions).
  This resolves #90.
- Changes how ellipsoids are computed/initialized so that all linear algebra
  is done before inverting the covariance matrix. This significantly improves
  eigenvector/eigenvalue behavior.
- The transformation axes are now the cholesky decomposition instead of the
  principal eigenvectors (and eigenvalues), although those are still stored and
  used at various stages. This leads to a moderate increase in overhead.
- Along with some small changes to how I've implemented `hslice`, this
  resolves #83.
- A number of calculations are now done in log-space rather than in linear
  space to avoid over/underflows. This resolves #91.

Slice sampling:
- Fixed sampling issues for `'rslice'` and `'hslice'`. The former now uses
  proposals from the unit n-sphere; the latter from the n-d standard normal.
  `'slice'` still uses the principle eigenvectors.

utils:
- Added `approx` option to `simulate_run` (should've been introduced before).

dynamicsampler:
- Fixed a bug where the sampler would return bounds (and related quantities)
  in `sampler.results` even when `save_bounds=False`. As before, bounds are
  still saved internally from `sample_initial` (since we need them to
  add batches).

demos:
- All the demos have been updated!

If all goes well, there shouldn't be any more large content updates
in the future other than #39.

.gitignore

@@ -48,3 +48,6 @@ docs/examples
 
 # old files
 old/*
+
+# saved data
+demos/*.pkl

docs/source/index.rst

@@ -12,28 +12,38 @@ Installation
 ============
 
 ``dynesty`` is compatible with both Python 2.7 and Python 3.6. It requires
-``numpy``, ``scipy``, ``matplotlib``, and ``six``. Installing the most
-recent stable version of the package is as easy as::
+``numpy``, ``scipy``, ``matplotlib``, and ``six``; ``sklearn`` is also
+recommended but not required. Installing the most recent stable version
+of the package is as easy as::
 
     pip install dynesty
 
 Alternately, for users who might want newer development versions, it can also
-be installed directly from the directory by running::
+be installed directly from a local copy of the repository by running::
 
     python setup.py install
 
 Changelog
 =========
 
-0.9.2 (2018-03-03)
+0.9.2 (2018-03-17)
 ------------------
 
-* Added in a fast approximation option for `jitter_run`.
+* Added in a fast approximation option for `jitter_run` and `simulate_run`.
 
-* Modified the default stopping heuristic. It now evaluates much faster but is
-  a less accurate probe of the "true" errors.
+* Modified the default stopping heuristic. It now evaluates significantly
+  faster but is a less accurate probe of the "true" KL divergence.
 
-* Added safety checks for `prior_transform`.
+* Modified `'rwalk'` behavior to better deal with edge cases.
+
+* Changed defaults so performance should now be more stable (albiet slower) 
+  for the average user.
+
+* Improved the stability of bounding ellipsoids.
+
+* Fixed performance issues with `'rslice'` and `'hslice'`.
+
+* Small plotting improvements.
 
 0.9.1 (2018-03-01)
 ------------------

dynesty/__init__.py

@@ -10,4 +10,4 @@
 from . import utils
 
 
-__version__ = "0.9.1e"
+__version__ = "0.9.2"

dynesty/bounding.py

@@ -33,25 +33,46 @@
 from scipy import special
 from scipy import spatial
 from scipy import linalg as lalg
-
-try:
-    from scipy.cluster.vq import kmeans2
-    HAVE_KMEANS = True
-except ImportError:  # pragma: no cover
-    HAVE_KMEANS = False
-
 from .utils import random_choice
 
 __all__ = ["UnitCube", "Ellipsoid", "MultiEllipsoid",
            "RadFriends", "SupFriends",
-           "vol_prefactor", "randsphere",
+           "vol_prefactor", "logvol_prefactor", "randsphere",
            "bounding_ellipsoid", "bounding_ellipsoids",
            "_bounding_ellipsoids", "_ellipsoid_bootstrap_expand",
            "_ellipsoids_bootstrap_expand", "_friends_bootstrap_radius",
            "_friends_leaveoneout_radius"]
 
 SQRTEPS = math.sqrt(float(np.finfo(np.float64).eps))
 
+# Try and import k-means clustering (used with 'multi').
+try:
+    from scipy.cluster.vq import kmeans2
+    HAVE_KMEANS = True
+except ImportError:  # pragma: no cover
+    HAVE_KMEANS = False
+
+# Alias for computing more robust/stable shrunk covariance estimates.
+try:
+    # Use the Oracle Approximating Shrinkage shrunk estimator.
+    # This is optimal for independent Normally-distributed data.
+    from sklearn.covariance import oas as oas_cov
+
+    def covariance(points):
+        return oas_cov(points)[0]
+    HAVE_OAS = True
+except ImportError:  # pragma: no cover
+    # If not available, use the default MLE sample covariance estimator.
+    # This is unbiased but noisy, especially if the covariance is sparse
+    # and the number of samples is small.
+    from numpy import cov as mle_cov
+
+    def covariance(points):
+        return mle_cov(points, rowvar=False)
+    HAVE_OAS = False
+    warnings.warn("Robust OAS shrunk covariance estimator from `sklearn` is "
+                  "not available. Defaulting to MLE estimator from `numpy`.")
+
 
 class UnitCube(object):
     """
@@ -137,36 +158,39 @@ class Ellipsoid(object):
     ctr : `~numpy.ndarray` with shape (N,)
         Coordinates of ellipsoid center.
 
-    am : `~numpy.ndarray` with shape (N, N)
-        Matrix describing the axes.
+    cov : `~numpy.ndarray` with shape (N, N)
+        Covariance matrix describing the axes.
 
     """
 
-    def __init__(self, ctr, am):
+    def __init__(self, ctr, cov):
         self.n = len(ctr)  # dimension
         self.ctr = np.array(ctr)  # center coordinates
-        self.am = np.array(am)  # precision matrix (inverse of covariance)
+        self.cov = np.array(cov)  # covariance matrix
+        self.am = lalg.pinvh(cov)  # precision matrix (inverse of covariance)
+        self.axes = lalg.cholesky(cov, lower=True)  # transformation axes
 
         # Volume of ellipsoid is the volume of an n-sphere divided
         # by the (determinant of the) Jacobian associated with the
         # transformation, which by definition is the precision matrix.
-        self.vol = vol_prefactor(self.n) / np.sqrt(linalg.det(self.am))
+        detsign, detln = linalg.slogdet(self.am)
+        self.vol = np.exp(logvol_prefactor(self.n) - 0.5 * detln)
 
         # The eigenvalues (l) of `a` are (a^-2, b^-2, ...) where
         # (a, b, ...) are the lengths of principle axes.
         # The eigenvectors (v) are the normalized principle axes.
-        l, v = linalg.eigh(self.am)
+        l, v = lalg.eigh(self.cov)
         if np.all((l > 0.) & (np.isfinite(l))):
-            self.axlens = 1. / np.sqrt(l)
+            self.axlens = np.sqrt(l)
         else:
             raise ValueError("The input precision matrix defining the "
                              "ellipsoid {0} is apparently singular with "
-                             "l={1} and v={2}.".format(self.am, l, v))
+                             "l={1} and v={2}.".format(self.cov, l, v))
 
-        # Scaled eigenvectors are the axes, where `axes[:,i]` is the
-        # i-th axis.  Multiplying this matrix by a vector will transform a
+        # Scaled eigenvectors are the principle axes, where `paxes[:,i]` is the
+        # i-th axis. Multiplying this matrix by a vector will transform a
         # point in the unit n-sphere to a point in the ellipsoid.
-        self.axes = np.dot(v, np.diag(self.axlens))
+        self.paxes = np.dot(v, np.diag(self.axlens))
 
         # Amount by which volume was increased after initialization (i.e.
         # cumulative factor from `scale_to_vol`).
@@ -175,8 +199,9 @@ def __init__(self, ctr, am):
     def scale_to_vol(self, vol):
         """Scale ellipoid to a target volume."""
 
-        f = (vol / self.vol) ** (1.0 / self.n)  # linear factor
+        f = np.exp((np.log(vol) - np.log(self.vol)) / self.n)  # linear factor
         self.expand *= f
+        self.cov *= f**2
         self.am *= f**-2
         self.axlens *= f
         self.axes *= f
@@ -186,7 +211,7 @@ def major_axis_endpoints(self):
         """Return the endpoints of the major axis."""
 
         i = np.argmax(self.axlens)  # find the major axis
-        v = self.axes[:, i]  # vector from center to major axis endpoint
+        v = self.paxes[:, i]  # vector from center to major axis endpoint
 
         return self.ctr - v, self.ctr + v
 
@@ -296,10 +321,12 @@ def update(self, points, pointvol=0., rstate=None, bootstrap=0,
         ell = bounding_ellipsoid(points, pointvol=pointvol)
         self.n = ell.n
         self.ctr = ell.ctr
+        self.cov = ell.cov
         self.am = ell.am
         self.vol = ell.vol
         self.axlens = ell.axlens
         self.axes = ell.axes
+        self.paxes = ell.paxes
         self.expand = ell.expand
 
         # Use bootstrapping to determine the volume expansion factor.
@@ -344,34 +371,36 @@ class MultiEllipsoid(object):
         Collection of coordinates of ellipsoid centers. Used to initialize
         :class:`MultiEllipsoid` if :data:`ams` is also provided.
 
-    ams : `~numpy.ndarray` with shape (M, N, N), optional
+    covs : `~numpy.ndarray` with shape (M, N, N), optional
         Collection of matrices describing the axes of the ellipsoids. Used to
         initialize :class:`MultiEllipsoid` if :data:`ctrs` also provided.
 
     """
 
-    def __init__(self, ells=None, ctrs=None, ams=None):
+    def __init__(self, ells=None, ctrs=None, covs=None):
         if ells is not None:
             # Try to initialize quantities using provided `Ellipsoid` objects.
-            if (ctrs is None) and (ams is None):
+            if (ctrs is None) and (covs is None):
                 self.nells = len(ells)
                 self.ells = ells
                 self.ctrs = np.array([ell.ctr for ell in self.ells])
+                self.covs = np.array([ell.cov for ell in self.ells])
                 self.ams = np.array([ell.am for ell in self.ells])
             else:
                 raise ValueError("You cannot specific both `ells` and "
                                  "(`ctrs`, `ams`)!")
         else:
             # Try to initialize quantities using provided `ctrs` and `ams`.
-            if (ctrs is None) and (ams is None):
+            if (ctrs is None) and (covs is None):
                 raise ValueError("You must specify either `ells` or "
-                                 "(`ctrs`, `ams`).")
+                                 "(`ctrs`, `covs`).")
             else:
                 self.nells = len(ctrs)
                 self.ctrs = np.array(ctrs)
-                self.ams = np.array(ams)
-                self.ells = [Ellipsoid(ctrs[i], ams[i])
+                self.covs = np.array(covs)
+                self.ells = [Ellipsoid(ctrs[i], covs[i])
                              for i in range(self.nells)]
+                self.ams = np.array([ell.am for ell in self.ells])
 
         # Compute quantities.
         self.vols = np.array([ell.vol for ell in self.ells])
@@ -394,9 +423,6 @@ def scale_to_vols(self, vols):
     def major_axis_endpoints(self):
         """Return the endpoints of the major axis of each ellipsoid."""
 
-        i = np.argmax(self.axlens)  # find the major axis
-        v = self.axes[:, i]  # vector from center to major axis endpoint
-
         return np.array([ell.major_axis_endpoints() for ell in self.ells])
 
     def within(self, x, j=None):
@@ -590,6 +616,7 @@ def update(self, points, pointvol=0., vol_dec=0.5, vol_check=2.,
         self.nells = len(ells)
         self.ells = ells
         self.ctrs = np.array([ell.ctr for ell in self.ells])
+        self.covs = np.array([ell.cov for ell in self.ells])
         self.ams = np.array([ell.am for ell in self.ells])
         self.vols = np.array([ell.vol for ell in self.ells])
         self.vol_tot = sum(self.vols)
@@ -666,7 +693,7 @@ def within(self, x, ctrs, kdtree=None):
             # If no K-D Tree is provided, execute a brute-force
             # search over all balls.
             nctrs = len(ctrs)
-            within = np.array([linalg.norm(ctrs[i] - x) <= self.radius
+            within = np.array([lalg.norm(ctrs[i] - x) <= self.radius
                                for i in range(nctrs)], dtype='bool')
             idxs = np.arange(nctrs)[within]
         else:
@@ -1116,14 +1143,32 @@ def vol_prefactor(n, p=2.):
     return f
 
 
+def logvol_prefactor(n, p=2.):
+    """
+    Returns the ln(volume constant) for an `n`-dimensional sphere with an
+    :math:`L^p` norm. The constant is defined as::
+
+        lnf = n * ln(2.) + n * LogGamma(1./p + 1) - LogGamma(n/p + 1.)
+
+    By default the `p=2.` norm is used (i.e. the standard Euclidean norm).
+
+    """
+
+    p *= 1.  # convert to float in case user inputs an integer
+    lnf = (n * np.log(2.) + n * special.gammaln(1./p + 1.) -
+           special.gammaln(n/p + 1))
+
+    return lnf
+
+
 def randsphere(n, rstate=None):
     """Draw a point uniformly within an `n`-dimensional unit sphere."""
 
     if rstate is None:
         rstate = np.random
 
     z = rstate.randn(n)  # initial n-dim vector
-    zhat = z / linalg.norm(z)  # normalize
+    zhat = z / lalg.norm(z)  # normalize
     xhat = zhat * rstate.rand()**(1./n)  # scale
 
     return xhat
@@ -1162,16 +1207,26 @@ def bounding_ellipsoid(points, pointvol=0.):
     if npoints == 1:
         if pointvol > 0.:
             ctr = points[0]
-            r = (pointvol / vol_prefactor(ndim))**(1./ndim)
-            am = (1. / r**2) * np.identity(ndim)
-            return Ellipsoid(ctr, am)
+            r = np.exp((np.log(pointvol) - logvol_prefactor(ndim)) / ndim)
+            covar = r**2 * np.identity(ndim)
+            return Ellipsoid(ctr, covar)
         else:
             raise ValueError("Cannot compute a bounding ellipsoid to a "
                              "single point if `pointvol` is not specified.")
 
     # Calculate covariance of points.
     ctr = np.mean(points, axis=0)
-    cov = np.cov(points, rowvar=False)
+    cov = covariance(points)
+
+    if not HAVE_OAS:
+        if npoints <= 2 * ndim:
+            warnings.warn("Volume is sparsely sampled. MLE covariance "
+                          "estimates and associated ellipsoid decompositions "
+                          "might be unstable.")
+        if ndim >= 50:
+            warnings.warn("Dimensionality is large. MLE covariance "
+                          "estimates and associated ellipsoid decompositions "
+                          "might be unstable.")
 
     # When ndim = 1, `np.cov` returns a 0-d array. Make it a 1x1 2-d array.
     if ndim == 1:
@@ -1189,20 +1244,20 @@ def bounding_ellipsoid(points, pointvol=0.):
     # nonsingular to deal with pathological cases where the ellipsoid has
     # "zero" volume. This can occur when `npoints <= ndim` or when enough
     # points are linear combinations of other points.
-    cov_temp = np.array(cov)
+    covar = np.array(cov)
     for trials in range(100):
         try:
             # Check if matrix is invertible.
-            am = lalg.solve(cov_temp, np.eye(ndim), assume_a='pos')
-            l, v = linalg.eigh(am)  # compute eigenvalues/vectors
+            am = lalg.pinvh(covar)
+            l, v = lalg.eigh(covar)  # compute eigenvalues/vectors
             if np.all((l > 0.) & (np.isfinite(l))):
                 break
             else:
                 raise RuntimeError("The eigenvalue/eigenvector decomposition "
                                    "failed!")
         except:
             # If the matrix remains singular/unstable, increase the nugget.
-            cov_temp = cov + np.eye(ndim) * (2.**(trials+1) * 1e-10)
+            covar = cov + np.eye(ndim) * (2.**(trials+1) * 1e-10)
             pass
     else:
         warnings.warn("Failed to guarantee the ellipsoid axes will be "
@@ -1223,16 +1278,16 @@ def bounding_ellipsoid(points, pointvol=0.):
     one_minus_a_bit = 1. - SQRTEPS
 
     if fmax > one_minus_a_bit:
-        am *= one_minus_a_bit / fmax
+        covar *= fmax / one_minus_a_bit
 
     # Initialize our ellipsoid.
-    ell = Ellipsoid(ctr, am)
+    ell = Ellipsoid(ctr, covar)
 
     # Expand our ellipsoid to encompass a minimum volume.
     if pointvol > 0.:
-        v = npoints * pointvol
-        if ell.vol < v:
-            ell.scale_to_vol(v)
+        minvol = npoints * pointvol
+        if ell.vol < minvol:
+            ell.scale_to_vol(minvol)
 
     return ell