Commit cdb30e2f authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Merge branch 'work-on-proposals' into 'master'

Testing improvements to the proposal

See merge request lscsoft/bilby!707
parents 9f4b7a4d 354cd4f3
......@@ -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):
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment