diff --git a/bilby/core/sampler/dynesty.py b/bilby/core/sampler/dynesty.py index b1d12fad48d379e6540497e233090cab02f3395f..7a4e49e4b931d98dcafeefd3dd2453aecca7150e 100644 --- a/bilby/core/sampler/dynesty.py +++ b/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):