A few typos and bugfixes

bilby/core/sampler/dynesty_utils.py

@@ -153,6 +153,7 @@ class BilbyMultiEllipsoidSampler(MultiEllipsoidSampler):
         else:
             i = self.rstate.randint(self.nlive)
             u = self.live_u[i, :]
+            # this should probably be ncdim once we are able to use that
             return u, np.identity(self.npdim)
 
     def update(self, pointvol):
@@ -266,6 +267,7 @@ def sample_rwalk_bilby(args):
 
     # Setup.
     n = len(u)
+    n_cluster = axes.shape[0]
     walks = kwargs.get("walks", 100)  # minimum number of steps
 
     axis_proposals = {"volumetric", "normal", "chi", "axis"}
@@ -282,7 +284,8 @@ def sample_rwalk_bilby(args):
             failures[prop] = 0
 
     if "live" in kwargs:
-        live = np.unique(deepcopy(kwargs["live"]), axis=0)
+        live = np.unique(deepcopy(kwargs["live"]), axis=0)[:, :n_cluster]
+        live = np.unique(live, axis=0)
         matches = np.where(np.equal(u, live).all(axis=1))[0]
         np.delete(live, matches, 0)
     else:
@@ -304,12 +307,12 @@ def sample_rwalk_bilby(args):
 
     for ii in range(walks):
         prop = props[ii % len(props)]
-        proposal_kwargs["u"] = u
-        proposed = _proposal_map[prop](**proposal_kwargs)
+        u_prop = np.random.uniform(0, 1, len(u))
+        proposed = _proposal_map[prop](u=u[:n_cluster], **proposal_kwargs)
         if isinstance(proposed, tuple):
-            u_prop, ln_jacobian = proposed
+            u_prop[:n_cluster], ln_jacobian = proposed
         else:
-            u_prop = proposed
+            u_prop[:n_cluster] = proposed
             ln_jacobian = 0
 
         if ln_jacobian < np.log(rstate.uniform(0, 1)):
@@ -317,12 +320,20 @@ def sample_rwalk_bilby(args):
             reject += 1
             continue
 
-        u_prop, fail = apply_boundary(
-            u_prop=u_prop,
-            periodic=periodic,
-            reflective=reflective,
-            nonbounded=nonbounded,
-        )
+        # Only apply boundary conditions if there is no jacobian for the
+        # proposal. This may be overly conservative, but some proposals
+        # that require a Jacobian, e.g., the stretch, break detailed
+        # balance if using a periodic or reflective boundary.
+        if ln_jacobian == 0:
+            u_prop, fail = apply_boundary(
+                u_prop=u_prop,
+                periodic=periodic,
+                reflective=reflective,
+                nonbounded=nonbounded
+            )
+        else:
+            fail = (np.min(u_prop) < 0) or (np.max(u_prop) > 1)
+
         if fail:
             failures[prop] += 1
             reject += 1
@@ -384,8 +395,8 @@ def propose_diff_evo(u, live, rstate, **kwargs):
     first, second = rstate.choice(nlive, 2, replace=False)
     diff = live[second] - live[first]
     if rstate.uniform(0, 1) < 0.5:
-        diff *= 2.38 / n ** 0.5
-        diff *= (100 ** rstate.uniform(0, 1)) / 10
+        diff *= 2.38 / n**0.5
+        diff *= (100**rstate.uniform(0, 1)) / 10
     u_prop = u + diff
     return u_prop
 
@@ -457,13 +468,13 @@ def propose_ensemble_stretch(u, live, rstate, **kwargs):
     """
     nlive, n = live.shape
     max_scale = 3
-    scale = np.exp(np.log(max_scale) * rstate.uniform(0, 2)) / max_scale
+    scale = ((max_scale - 1.0) * rstate.uniform() + 1) ** 2.0 / max_scale
 
     other = rstate.choice(nlive)
     other = live[other]
 
     u_prop = other + scale * (u - other)
-    ln_jacobian = np.log(scale) * n
+    ln_jacobian = np.log(scale) * (n - 1)
     return u_prop, ln_jacobian