Repository: joshspeagle/dynesty
Issue: 1.0.2 issue with rejecting invalid multiellipsoid region

Hi,

while running the latest master version (with the 'rejecting multiellipsoid region' change), I'm getting this exception. I didn't have it before (with the previously released version). 
The error seems to occur exactly after the warnign about rejecting the multiellipsoid region.

```
eff(%):  0.282 | loglstar: 10019.867 < 10026.423 < 10028.004 | logz: 9949.172 +/-  0.545 | stop:  1.314]
/home/skoposov/pyenv38/lib/python3.8/site-packages/dynesty/bounding.py:596: UserWarning: Rejecting invalid MultiEllipsoid region
  warnings.warn('Rejecting invalid MultiEllipsoid region')
59036it [8:18:40,  1.97it/s, batch: 8 | bound: 787 | nc: 332 | ncall: 20754337 | eff(%):  0.282 | loglstar: 10019.867 < 10026.430 < 10028.004 | logz: 9949.172 +/-  0.545 | stop:  1.314]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-2-4fd4a960409f> in <module>
----> 1 R=qq.sampler()

~/science/dipper/qq.py in sampler()
    359                                             sample='slice'
    360                                             )
--> 361     dsampler.run_nested(dlogz_init=0.01,
    362                         nlive_init=500,
    363                         nlive_batch=500,

~/pyenv38/lib/python3.8/site-packages/dynesty/dynamicsampler.py in run_nested(self, nlive_init, maxiter_init, maxcall_init, dlogz_init, logl_max_init, n_effective_init, nlive_batch, wt_function, wt_kwargs, maxiter_batch, maxcall_batch, maxiter, maxcall, maxbatch, n_effective, stop_function, stop_kwargs, use_stop, save_bounds, print_progress, print_func, live_points)
   1658                     # Compute our sampling bounds using the provided
   1659                     # weight function.
-> 1660                     passback = self.add_batch(nlive=nlive_batch,
   1661                                               wt_function=wt_function,
   1662                                               wt_kwargs=wt_kwargs,

~/pyenv38/lib/python3.8/site-packages/dynesty/dynamicsampler.py in add_batch(self, nlive, wt_function, wt_kwargs, maxiter, maxcall, logl_bounds, save_bounds, print_progress, print_func, stop_val)
   1767                 if logl_bounds is None:
   1768                     logl_bounds = wt_function(res, wt_kwargs)
-> 1769                 for results in self.sample_batch(nlive_new=nlive,
   1770                                                  logl_bounds=logl_bounds,
   1771                                                  maxiter=maxiter,

~/pyenv38/lib/python3.8/site-packages/dynesty/dynamicsampler.py in sample_batch(self, nlive_new, update_interval, logl_bounds, maxiter, maxcall, save_bounds)
   1207         # sample past the original bounds "for free".
   1208         for i in range(1):
-> 1209             for it, results in enumerate(self.sampler.sample(dlogz=0.,
   1210                                          logl_max=logl_max,
   1211                                          maxiter=maxiter-nlive_new-1,

~/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py in sample(self, maxiter, maxcall, dlogz, logl_max, n_effective, add_live, save_bounds, save_samples)
    780             # `logl > loglstar` using the bounding distribution and sampling
    781             # method from our sampler.
--> 782             u, v, logl, nc = self._new_point(loglstar_new, logvol)
    783             ncall += nc
    784             self.ncall += nc

~/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py in _new_point(self, loglstar, logvol)
    378         while True:
    379             # Get the next point from the queue
--> 380             u, v, logl, nc, blob = self._get_point_value(loglstar)
    381             ncall += nc
    382 

~/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py in _get_point_value(self, loglstar)
    362         # If the queue is empty, refill it.
    363         if self.nqueue <= 0:
--> 364             self._fill_queue(loglstar)
    365 
    366         # Grab the earliest entry.

~/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py in _fill_queue(self, loglstar)
    331             if self._beyond_unit_bound(loglstar):
    332                 # Propose points using the provided sampling/bounding options.
--> 333                 point, axes = self.propose_point()
    334                 evolve_point = self.evolve_point
    335             else:

~/pyenv38/lib/python3.8/site-packages/dynesty/nestedsamplers.py in propose_live(self)
    664 
    665         # Pick a random ellipsoid that encompasses `u`.
--> 666         ell_idx = ell_idxs[self.rstate.randint(nidx)]
    667 
    668         # Choose axes.

mtrand.pyx in numpy.random.mtrand.RandomState.randint()

_bounded_integers.pyx in numpy.random._bounded_integers._rand_int64()

ValueError: low >= high
```
I'm not sure I'll be able to provide a backtrace on this one, as it  was a very long run.


Changed files:
dynesty/bounding.py

diff --git a/dynesty/bounding.py b/dynesty/bounding.py
index 140855c..36223b1 100644
--- a/dynesty/bounding.py
+++ b/dynesty/bounding.py
@@ -593,8 +593,7 @@ class MultiEllipsoid(object):
         # Sanity check: all points must be contained in some ellipsoid
         if not all(any(ell.contains(p) for ell in ells) for p in points):
             # refuse to update
-            warnings.warn('Rejecting invalid MultiEllipsoid region')
-            return
+            raise RuntimeError('Rejecting invalid MultiEllipsoid region')
 
         # Update the set of ellipsoids.
         self.nells = len(ells)
@@ -1364,9 +1363,13 @@ def bounding_ellipsoid(points, pointvol=0.):
     f = np.einsum('...i, ...i', np.tensordot(delta, am, axes=1), delta)
     fmax = np.max(f)
 
+    
     # Due to round-off errors, we actually scale the ellipsoid so the
     # outermost point obeys `(x-v)^T A (x-v) < 1 - (a bit) < 1`.
-    one_minus_a_bit = 1. - SQRTEPS
+    ROUND_DELTA = 1e-5
+    # numerical experiments show that round off errors can reach large
+    # values if the matrix eigen values are very low
+    one_minus_a_bit = 1. - ROUND_DELTA
     covar_mod = np.array(covar)
     if fmax > one_minus_a_bit:
         covar_mod *= fmax / one_minus_a_bit
@@ -1375,6 +1378,7 @@ def bounding_ellipsoid(points, pointvol=0.):
     # us numerically unstable again. Just to be **ultra safe**.
     covar2 = np.array(covar_mod)
     for trials in range(ntries):
+        failed = False
         try:
             # Check if matrix is invertible.
             am = lalg.pinvh(covar2)
@@ -1387,19 +1391,23 @@ def bounding_ellipsoid(points, pointvol=0.):
             if np.all((l > 0.) & (np.isfinite(l))):
                 break
             else:
-                raise RuntimeError("The eigenvalue/eigenvector decomposition "
-                                   "failed!")
-        except:
+                failed =True
+        except lalg.LinAlgError:
+            failed =True
+        if failed:
             # If the matrix remains singular/unstable,
             # suppress the off-diagonal elements.
             coeff = 1.1**(trials+1) / 1.1**100
             covar2 = (1. - coeff) * covar_mod + coeff * np.eye(ndim)
-            pass
     else:
         warnings.warn("Failed to guarantee the ellipsoid axes will be "
                       "non-singular. Defaulting to last working axes.")
         covar2 = np.array(covar)  # default to last safe version
-
+    # this is a final check
+    # if this fails the ellipsoid is broken already
+    if np.einsum('...i, ...i', np.tensordot(delta, lalg.pinvh(covar2), axes=1), delta).max()>=1:
+        raise RuntimeError("Failed to initialize the ellipsoid to contain all the points")
+        
     # Initialize our ellipsoid with *safe* covariance matrix.
     ell = Ellipsoid(ctr, covar2)
 
@@ -1458,54 +1466,51 @@ def _bounding_ellipsoids(points, ell, pointvol=0., vol_dec=0.5,
     start_ctrs = np.vstack((p1, p2))  # shape is (k, ndim) = (2, ndim)
 
     # Split points into two clusters using k-means clustering with k=2.
-    try:
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore")
-            k2_res = kmeans2(points, k=start_ctrs, iter=10, minit='matrix',
-                             check_finite=False)
-        labels = k2_res[1]  # cluster identifier ; shape is (npoints,)
-
-        # Get points in each cluster.
-        points_k = [points[labels == k, :] for k in (0, 1)]
-
-        # If either cluster has less than ndim+1 points, the bounding ellipsoid
-        # will be ill-constrained. Reject the split and simply return the
-        # original ellipsoid bounding all the points.
-        if points_k[0].shape[0] < 2 * ndim or points_k[1].shape[0] < 2 * ndim:
-            return [ell]
-
-        # Bounding ellipsoid for each cluster, possibly enlarged
-        # to a minimum volume.
-        ells = [bounding_ellipsoid(points_j, pointvol=pointvol)
-                for points_j in points_k]
-
-        # If the total volume decreased by a factor of `vol_dec`, we accept
-        # the split into subsets. We then recursively split each subset.
-        if ells[0].vol + ells[1].vol < vol_dec * ell.vol:
-            return (_bounding_ellipsoids(points_k[0], ells[0],
-                                         pointvol=pointvol, vol_dec=vol_dec,
-                                         vol_check=vol_check) +
-                    _bounding_ellipsoids(points_k[1], ells[1],
-                                         pointvol=pointvol, vol_dec=vol_dec,
-                                         vol_check=vol_check))
-
-        # Otherwise, see if the total ellipsoid volume is larger than the
-        # minimum volume by a factor of `vol_check`. If it is, this indicates
-        # that there may be more than 2 clusters and we should try to
-        # subdivide further.
-        if ell.vol > vol_check * npoints * pointvol:
-            out = (_bounding_ellipsoids(points_k[0], ells[0],
-                                        pointvol=pointvol, vol_dec=vol_dec,
-                                        vol_check=vol_check) +
-                   _bounding_ellipsoids(points_k[1], ells[1],
-                                        pointvol=pointvol, vol_dec=vol_dec,
-                                        vol_check=vol_check))
-
-            # Only accept the split if the volume decreased significantly.
-            if sum(e.vol for e in out) < vol_dec * ell.vol:
-                return out
-    except:
-        pass
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        k2_res = kmeans2(points, k=start_ctrs, iter=10, minit='matrix',
+                         check_finite=False)
+    labels = k2_res[1]  # cluster identifier ; shape is (npoints,)
+
+    # Get points in each cluster.
+    points_k = [points[labels == k, :] for k in (0, 1)]
+
+    # If either cluster has less than ndim+1 points, the bounding ellipsoid
+    # will be ill-constrained. Reject the split and simply return the
+    # original ellipsoid bounding all the points.
+    if points_k[0].shape[0] < 2 * ndim or points_k[1].shape[0] < 2 * ndim:
+        return [ell]
+
+    # Bounding ellipsoid for each cluster, possibly enlarged
+    # to a minimum volume.
+    ells = [bounding_ellipsoid(points_j, pointvol=pointvol)
+            for points_j in points_k]
+
+    # If the total volume decreased by a factor of `vol_dec`, we accept
+    # the split into subsets. We then recursively split each subset.
+    if ells[0].vol + ells[1].vol < vol_dec * ell.vol:
+        return (_bounding_ellipsoids(points_k[0], ells[0],
+                                     pointvol=pointvol, vol_dec=vol_dec,
+                                     vol_check=vol_check) +
+                _bounding_ellipsoids(points_k[1], ells[1],
+                                     pointvol=pointvol, vol_dec=vol_dec,
+                                     vol_check=vol_check))
+
+    # Otherwise, see if the total ellipsoid volume is larger than the
+    # minimum volume by a factor of `vol_check`. If it is, this indicates
+    # that there may be more than 2 clusters and we should try to
+    # subdivide further.
+    if ell.vol > vol_check * npoints * pointvol:
+        out = (_bounding_ellipsoids(points_k[0], ells[0],
+                                    pointvol=pointvol, vol_dec=vol_dec,
+                                    vol_check=vol_check) +
+               _bounding_ellipsoids(points_k[1], ells[1],
+                                    pointvol=pointvol, vol_dec=vol_dec,
+                                    vol_check=vol_check))
+
+        # Only accept the split if the volume decreased significantly.
+        if sum(e.vol for e in out) < vol_dec * ell.vol:
+            return out
 
     # Otherwise, we are happy with the single bounding ellipsoid.
     return [ell]