Testing improvements to the proposal

bilby/core/sampler/dynesty.py

@@ -17,7 +17,6 @@ from .base_sampler import Sampler, NestedSampler
 from numpy import linalg
 from dynesty.utils import unitcheck
 import warnings
-import math
 
 
 class Dynesty(NestedSampler):
@@ -559,6 +558,7 @@ def sample_rwalk_bilby(args):
     walks = kwargs.get('walks', 25)  # minimum number of steps
     maxmcmc = kwargs.get('maxmcmc', 2000)  # Maximum number of steps
     nact = kwargs.get('nact', 5)  # Number of ACT
+    old_act = kwargs.get('old_act', walks)
 
     # Initialize internal variables
     accept = 0
@@ -573,28 +573,12 @@ def sample_rwalk_bilby(args):
     drhat, dr, du, u_prop, logl_prop = np.nan, np.nan, np.nan, np.nan, np.nan
     while len(u_list) < nact * act:
 
-        if scale == 0.:
-            raise RuntimeError("The random walk sampling is stuck! "
-                               "Some useful output quantities:\n"
-                               "u: {0}\n"
-                               "drhat: {1}\n"
-                               "dr: {2}\n"
-                               "du: {3}\n"
-                               "u_prop: {4}\n"
-                               "loglstar: {5}\n"
-                               "logl_prop: {6}\n"
-                               "axes: {7}\n"
-                               "scale: {8}."
-                               .format(u, drhat, dr, du, u_prop,
-                                       loglstar, logl_prop, axes, scale))
-
         # Propose a direction on the unit n-sphere.
         drhat = rstate.randn(n)
         drhat /= linalg.norm(drhat)
 
         # Scale based on dimensionality.
-        # dr = drhat * rstate.rand()**(1. / n)  # CHANGED FROM DYNESTY 1.0
-        dr = drhat * rstate.rand(n)
+        dr = drhat * rstate.rand()**(1. / n)
 
         # Transform to proposal distribution.
         du = np.dot(axes, dr)
@@ -619,14 +603,6 @@ def sample_rwalk_bilby(args):
                 logl_list.append(logl_list[-1])
             continue
 
-        # Check if we're stuck generating bad numbers.
-        if nfail > 100 * walks:
-            warnings.warn("Random number generation appears to be "
-                          "extremely inefficient. Adjusting the "
-                          "scale-factor accordingly.")
-            nfail = 0
-            scale *= math.exp(-1. / n)
-
         # Check proposed point.
         v_prop = prior_transform(np.array(u_prop))
         logl_prop = loglikelihood(np.array(v_prop))
@@ -649,7 +625,8 @@ def sample_rwalk_bilby(args):
         # If we've taken the minimum number of steps, calculate the ACT
         if accept + reject > walks:
             act = estimate_nmcmc(
-                accept_ratio=accept / (accept + reject + nfail), maxmcmc=maxmcmc)
+                accept_ratio=accept / (accept + reject + nfail),
+                old_act=old_act, maxmcmc=maxmcmc)
 
         # If we've taken too many likelihood evaluations then break
         if accept + reject > maxmcmc and accept > 0:
@@ -663,16 +640,9 @@ def sample_rwalk_bilby(args):
                 # Break if we are above maxmcmc and have at least one accepted point
                 break
 
-        # Check if we're stuck generating bad points.
-        if accept + reject > 50 * walks:
-            scale *= math.exp(-1. / n)
-            warnings.warn("Random walk proposals appear to be "
-                          "extremely inefficient. Adjusting the "
-                          "scale-factor accordingly.")
-
     # If the act is finite, pick randomly from within the chain
-    if np.isfinite(act) and act < len(u_list):
-        idx = np.random.randint(act, len(u_list))
+    if np.isfinite(act) and int(.5 * nact * act) < len(u_list):
+        idx = np.random.randint(int(.5 * nact * act), len(u_list))
         u = u_list[idx]
         v = v_list[idx]
         logl = logl_list[idx]
@@ -689,12 +659,13 @@ def sample_rwalk_bilby(args):
         logl = logl_list[idx]
 
     blob = {'accept': accept, 'reject': reject, 'fail': nfail, 'scale': scale}
+    kwargs["old_act"] = act
 
     ncall = accept + reject
     return u, v, logl, ncall, blob
 
 
-def estimate_nmcmc(accept_ratio, maxmcmc, safety=5, tau=None):
+def estimate_nmcmc(accept_ratio, old_act, maxmcmc, safety=5, tau=None):
     """ Estimate autocorrelation length of chain using acceptance fraction
 
     Using ACL = (2/acc) - 1 multiplied by a safety margin. Code adapated from
@@ -706,8 +677,8 @@ def estimate_nmcmc(accept_ratio, maxmcmc, safety=5, tau=None):
     ----------
     accept_ratio: float [0, 1]
         Ratio of the number of accepted points to the total number of points
-    minmcmc: int
-        The minimum length of the MCMC chain to use
+    old_act: int
+        The ACT of the last iteration
     maxmcmc: int
         The maximum length of the MCMC chain to use
     safety: int
@@ -720,12 +691,14 @@ def estimate_nmcmc(accept_ratio, maxmcmc, safety=5, tau=None):
         tau = maxmcmc / safety
 
     if accept_ratio == 0.0:
-        return np.inf
+        Nmcmc_exact = (1 + 1 / tau) * old_act
     else:
-        Nmcmc_exact = (safety / tau) * (2. / accept_ratio - 1.)
+        Nmcmc_exact = (
+            (1. - 1. / tau) * old_act +
+            (safety / tau) * (2. / accept_ratio - 1.)
+        )
         Nmcmc_exact = float(min(Nmcmc_exact, maxmcmc))
-        Nmcmc = max(safety, int(Nmcmc_exact))
-        return Nmcmc
+    return max(safety, int(Nmcmc_exact))
 
 
 class DynestySetupError(Exception):