remove the feature of reverting the step sampler.

Old behaviour:
When running with MPI, step samplers can progress at different efficiencies.
The point of the first one that finishes N steps is accepted and this increases
the nested sampling likelihood threshold. The other step samplers are informed
and may revert / back track to an earlier state. This can cause some areas of
the parameter space to be explored much slower, and thus become undersampled, leading to biases.
see #159

New behaviour:
Reverting is disabled. The step sampler remembers the Lmin it was started
with, and returns the Lmin->L edge to the ReactiveNestedSampling integrator,
which stores it and may use it (unless the new Lmin is already above L).

This should reduce the issue, but is still not perfect, as the edge may
arrive too late. Ideally, we should wait until all stepsamplers at the
current Lmin have finished, but this is difficult to synchronize.

Resuming could pick up the late edges; however, we would need to sort
the point store by Lmin. This should be technically possible within
an improvement loop if we keep track of the last pointstore entry
at the start of the loop, and keep the entries since then sorted.

setup.cfg

@@ -2,7 +2,7 @@
 style = numpy
 check-return-types = False
 exclude = docs
-extend-ignore = E501,F401,E128,E231,E124,SIM114,DOC105,DOC106,DOC107,DOC301,DOC501,DOC503,DOC203,B006,SIM102,SIM113,DOC202
+extend-ignore = E501,F401,E128,E231,E124,SIM114,DOC105,DOC106,DOC107,DOC301,DOC501,DOC503,DOC203,B006,SIM102,SIM113,DOC202,DOC403,DOC404
 per-file-ignores =
     ultranest/plot.py: B006
     ultranest/integrator.py: B006

tests/test_popstepsampling.py

@@ -245,7 +245,7 @@ def test_SimpleSliceSampler_SLOW(seed=4):
     # resetting the seed to check the slice axes
     np.random.seed(seed)
     for i in range(popsize):
-        u[i],_,L[i],_= stepsampler.__next__(region, Lmin, us.copy(), Ls.copy(), transform, loglike_vectorized, test=True)
+        u[i],_,L[i],_,_= stepsampler.__next__(region, Lmin, us.copy(), Ls.copy(), transform, loglike_vectorized, test=True)
 
     # Basic check
     assert (L>Lmin).all(), (L,Lmin) # Lmin check

tests/test_stepsampling.py

@@ -224,16 +224,18 @@ def test_stepsampler(plot=False):
 
     stepsampler = CubeMHSampler(nsteps=len(paramnames))
     while True:
-        u1, p1, L1, nc = stepsampler.__next__(region, -1e100, region.u, Ls, transform, loglike)
+        u1, p1, L1, Lmin, nc = stepsampler.__next__(region, -1e100, region.u, Ls, transform, loglike)
         if u1 is not None:
             break
     assert L1 > -1e100
+    assert Lmin == -1e100
     print(u1, L1)
     while True:
-        u2, p2, L2, nc = stepsampler.__next__(region, -1e100, region.u, Ls, transform, loglike)
+        u2, p2, L2, Lmin, nc = stepsampler.__next__(region, -1e100, region.u, Ls, transform, loglike)
         if u2 is not None:
             break
     assert L2 > -1e100
+    assert Lmin == -1e100
     print(u2, L2)
     assert np.all(u1 != u2)
     assert np.all(L1 != L2)
@@ -253,7 +255,7 @@ def test_stepsampler_adapt_when_stuck(plot=False):
     for i in range(1000):
         if i > 100:
             assert False, i
-        unew, pnew, Lnew, nc = stepsampler.__next__(region, Lmin, us, Ls, transform, loglike, ndraw=10)
+        unew, pnew, Lnew, Lnew, nc = stepsampler.__next__(region, Lmin, us, Ls, transform, loglike, ndraw=10)
         if unew is not None:
             break
 
@@ -269,7 +271,7 @@ def test_stepsampler_adapt_when_stuck(plot=False):
         for i in range(1000):
             if i > 100:
                 assert False, i
-            unew, pnew, Lnew, nc = stepsampler.__next__(region, Lmin, us, Ls, transform, loglike, ndraw=10)
+            unew, pnew, Lnew, Lnew, nc = stepsampler.__next__(region, Lmin, us, Ls, transform, loglike, ndraw=10)
             if unew is not None:
                 break
 
@@ -301,7 +303,7 @@ def test_stepsampler_adapt(plot=True):
                 old_scale = stepsampler.scale
                 for i in range(5):
                     while True:
-                        unew, pnew, Lnew, nc = stepsampler.__next__(region, -1e100, region.u, Ls, transform, loglike)
+                        unew, pnew, Lnew, Lnew, nc = stepsampler.__next__(region, -1e100, region.u, Ls, transform, loglike)
                         if unew is not None:
                             break
                 new_scale = stepsampler.scale

ultranest/calibrator.py

@@ -80,17 +80,13 @@ def __init__(self,
 
         Parameters
         ----------
-
         param_names: list of str
             Names of the parameters.
             Length gives dimensionality of the sampling problem.
-
         loglike: function
             log-likelihood function.
-
         transform: function
             parameter transform from unit cube to physical parameters.
-
         kwargs: dict
             further arguments passed to ReactiveNestedSampler.
             if `log_dir` is set, then the suffix `-nsteps%d` is added for each

ultranest/integrator.py

@@ -14,6 +14,7 @@
 from __future__ import division, print_function
 
 import csv
+import itertools
 import json
 import operator
 import os
@@ -1893,7 +1894,7 @@ def _create_point(self, Lmin, ndraw, active_u, active_values):
             while ib >= len(self.samples):
                 ib = 0
                 if use_stepsampler:
-                    u, v, logl, nc = self.stepsampler.__next__(
+                    u, v, logl, Lmin_sampled, nc = self.stepsampler.__next__(
                         self.region,
                         transform=self.transform, loglike=self.loglike,
                         Lmin=Lmin, us=active_u, Ls=active_values,
@@ -1926,16 +1927,22 @@ def _create_point(self, Lmin, ndraw, active_u, active_values):
                     self.samplesv = np.concatenate(recv_samplesv, axis=0)
                     self.likes = np.concatenate(recv_likes, axis=0)
                     self.ncall += sum(recv_nc)
+                    if use_stepsampler:
+                        recv_Lmin_sampled = self.comm.gather(Lmin_sampled, root=0)
+                        self.Lmin_sampled = self.comm.bcast(recv_Lmin_sampled, root=0)
+                    else:
+                        self.Lmin_sampled = itertools.repeat(Lmin)
                 else:
                     self.samples = u
                     self.samplesv = v
                     self.likes = logl
                     self.ncall += nc
+                    self.Lmin_sampled = itertools.repeat(Lmin)
 
                 if self.log:
-                    for ui, vi, logli in zip(self.samples, self.samplesv, self.likes):
+                    for ui, vi, logli, Lmini in zip(self.samples, self.samplesv, self.likes, self.Lmin_sampled):
                         self.pointstore.add(
-                            _listify([Lmin, logli, quality], ui, vi),
+                            _listify([Lmini, logli, quality], ui, vi),
                             self.ncall)
 
             if self.likes[ib] > Lmin:

ultranest/mlfriends.pyx

@@ -642,7 +642,7 @@ class AffineLayer(ScalingLayer):
         ctr: vector
             Center of points
         T: matrix
-            Transformation matrix. This matrix whitens the points 
+            Transformation matrix. This matrix whitens the points
             to a unit Gaussian.
         invT: matrix
             Inverse transformation matrix. For transforming a unit
@@ -701,7 +701,7 @@ class AffineLayer(ScalingLayer):
         # Transformation matrix with the correct scale
         # this matrix whitens the points to a unit Gaussian.
         self.T = eigvec * eigval**-0.5
-        # Inverse transformation matrix, for transforming a unit 
+        # Inverse transformation matrix, for transforming a unit
         # Gaussian into something with the sample cov.
         self.invT = np.linalg.inv(self.T)
         # These also are the principle axes of the space
@@ -755,15 +755,15 @@ class MaxPrincipleGapAffineLayer(AffineLayer):
     """Affine whitening transformation.
 
     For learning the next layer's covariance, the clustering
-    and principal axis is considered: 
+    and principal axis is considered:
     the sample covariance is computed after subtracting
     the cluster mean. All points are projected onto the line
     defined by the principle axis vector starting from the origin.
     Then, on the sorted line positions, the largest gap is identified.
     All points before the gap are mean-subtracted, and all points
     after the gap are mean-subtracted. Then, the final
     sample covariance is computed. This should give a more "local"
-    covariance, even in the case where clusters could not yet be 
+    covariance, even in the case where clusters could not yet be
     clearly identified.
     """
 
@@ -819,7 +819,7 @@ class MaxPrincipleGapAffineLayer(AffineLayer):
 class LocalAffineLayer(AffineLayer):
     """Affine whitening transformation.
 
-    For learning the next layer's covariance, the points within 
+    For learning the next layer's covariance, the points within
     the MLradius are co-centered. This should give a more "local"
     covariance.
     """

ultranest/popstepsampler.py

@@ -355,7 +355,7 @@ def __next__(
             nc = 0
 
         u, p, L = self.prepared_samples.pop(0)
-        return u, p, L, nc
+        return u, p, L, Lmin, nc
 
 
 class PopulationSliceSampler(GenericPopulationSampler):
@@ -395,6 +395,7 @@ def __init__(
         self.scale_adapt_factor = scale_adapt_factor
         self.allu = []
         self.allL = []
+        self.allLmin = []
         self.currentt = []
         self.currentv = []
         self.currentp = []
@@ -428,6 +429,7 @@ def _setup(self, ndim):
         """Allocate arrays."""
         self.allu = np.zeros((self.popsize, self.nsteps + 1, ndim)) + np.nan
         self.allL = np.zeros((self.popsize, self.nsteps + 1)) + np.nan
+        self.allLmin = np.zeros((self.popsize, self.nsteps + 1)) + np.nan
         self.currentt = np.zeros(self.popsize) + np.nan
         self.currentv = np.zeros((self.popsize, ndim)) + np.nan
         self.generation = np.zeros(self.popsize, dtype=int_dtype) - 1
@@ -436,7 +438,7 @@ def _setup(self, ndim):
         self.searching_left = np.zeros(self.popsize, dtype=bool)
         self.searching_right = np.zeros(self.popsize, dtype=bool)
 
-    def setup_start(self, us, Ls, starting):
+    def setup_start(self, us, Ls, Lmin, starting):
         """Initialize walker starting points.
 
         For iteration zero, randomly selects a live point as starting point.
@@ -447,6 +449,8 @@ def setup_start(self, us, Ls, starting):
             live points
         Ls: np.array(nlive)
             loglikelihoods live points
+        Lmin: float
+            loglikelihood threshold
         starting: np.array(nwalkers, dtype=bool)
             which walkers to initialize.
 
@@ -466,6 +470,7 @@ def setup_start(self, us, Ls, starting):
 
         self.allu[starting,0] = us[i]
         self.allL[starting,0] = Ls[i]
+        self.allLmin[starting,:] = Lmin
         self.generation[starting] = 0
 
     @property
@@ -506,7 +511,7 @@ def _setup_currentp(self, nparams):
             print("setting currentp")
         self.currentp = np.zeros((self.popsize, nparams)) + np.nan
 
-    def advance(self, transform, loglike, Lmin, region):
+    def advance(self, transform, loglike, region):
         """Advance the walker population.
 
         Parameters
@@ -515,8 +520,6 @@ def advance(self, transform, loglike, Lmin, region):
             prior transform function
         loglike: function
             loglikelihood function
-        Lmin: float
-            current log-likelihood threshold
         region: MLFriends object
             Region
 
@@ -533,6 +536,7 @@ def advance(self, transform, loglike, Lmin, region):
             args = [
                 self.allu[i, self.generation],
                 self.allL[i, self.generation],
+                self.allLmin[i, self.generation],
                 # pass values directly
                 self.currentt,
                 self.currentv,
@@ -546,6 +550,7 @@ def advance(self, transform, loglike, Lmin, region):
             args = [
                 self.allu[movable, self.generation[movable]],
                 self.allL[movable, self.generation[movable]],
+                self.allLmin[movable, self.generation[movable]],
                 # this makes copies
                 self.currentt[movable],
                 self.currentv[movable],
@@ -556,16 +561,18 @@ def advance(self, transform, loglike, Lmin, region):
             ]
         if self.log:
             print("evolve will advance:", movable)
-
         uorig = args[0].copy()
         (
             (
                 currentt, currentv,
                 current_left, current_right, searching_left, searching_right
             ),
-            (success, unew, pnew, Lnew),
+            (success, unew, pnew, Lnew, Lminnew),
             nc
-        ) = evolve(transform, loglike, Lmin, *args)
+        ) = evolve(transform, loglike, *args)
+
+        if self.log:
+            print("evolve moved:", args[1], args[2], Lnew, Lminnew)
 
         if success.any():
             far_enough, (move_distance, reference_distance) = diagnose_move_distances(region, uorig[success,:], unew)
@@ -661,12 +668,14 @@ def __next__(
         if len(self.allu) == 0:
             self._setup(ndim)
 
-        step_back(Lmin, self.allL, self.generation, self.currentt)
+        # step_back(Lmin, self.allL, self.generation, self.currentt)
 
         starting = self.generation < 0
+        assert np.isfinite(Lmin), Lmin
         if starting.any():
-            self.setup_start(us[Ls > Lmin], Ls[Ls > Lmin], starting)
+            self.setup_start(us[Ls > Lmin], Ls[Ls > Lmin], Lmin, starting)
         assert (self.generation >= 0).all(), self.generation
+        assert np.isfinite(self.allLmin).all(), self.allLmin
 
         # find those where bracket is undefined:
         mask_starting = ~np.isfinite(self.currentt)
@@ -675,30 +684,35 @@ def __next__(
 
         if self.log:
             print(str(self), "(before)")
-        nc = self.advance(transform, loglike, Lmin, region)
+        nc = self.advance(transform, loglike, region)
         if self.log:
             print(str(self), "(after)")
+        assert np.isfinite(self.allLmin).all(), self.allLmin
 
         # harvest top individual if possible
         if self.generation[self.ringindex] == self.nsteps:
             if self.log:
                 print("have a candidate")
-            u, p, L = self.allu[self.ringindex, self.nsteps, :].copy(), self.currentp[self.ringindex, :].copy(), self.allL[self.ringindex, self.nsteps].copy()
+            u = self.allu[self.ringindex, self.nsteps, :].copy()
+            p = self.currentp[self.ringindex, :].copy()
+            L = self.allL[self.ringindex, self.nsteps].copy()
+            Lmin_start = self.allLmin[self.ringindex, self.nsteps].copy()
             assert np.isfinite(u).all(), u
             assert np.isfinite(p).all(), p
             self.generation[self.ringindex] = -1
             self.currentt[self.ringindex] = np.nan
             self.allu[self.ringindex,:,:] = np.nan
             self.allL[self.ringindex,:] = np.nan
+            self.allLmin[self.ringindex,:] = np.nan
 
             # adjust guess length
             newscale = (self.current_right[self.ringindex] - self.current_left[self.ringindex]) / 2
             self.scale = self.scale * 0.9 + 0.1 * newscale
 
             self.shift()
-            return u, p, L, nc
+            return u, p, L, Lmin_start, nc
         else:
-            return None, None, None, nc
+            return None, None, None, None, nc
 
 
 def slice_limit_to_unitcube(tleft, tright):
@@ -997,7 +1011,7 @@ def __next__(
             nc = 0
 
         u, p, L = self.prepared_samples.pop(0)
-        return u, p, L, nc
+        return u, p, L, Lmin, nc
 
 
 __all__ = [