DEV: add distance ratio as part of distance insertion test

Author: ColmTalbot

py/dynesty/dynamicsampler.py

@@ -605,6 +605,7 @@ def _configure_batch_sampler(main_sampler,
     saved_scale = np.array(main_sampler.saved_run['scale'])
     saved_blobs = np.array(main_sampler.saved_run['blob'])
     saved_distances_indices = np.array(main_sampler.saved_run["distance_insertion_index"])
+    saved_log_distance_ratios = np.array(main_sampler.saved_run["log_distance_ratio"])
     saved_likelihood_indices = np.array(main_sampler.saved_run["likelihood_insertion_index"])
     first_points = []
 
@@ -724,6 +725,7 @@ def _configure_batch_sampler(main_sampler,
 
         live_scale = saved_scale[subset0[0]]
         live_distance_index = saved_distances_indices[subset0[0]]
+        live_log_distance_ratio = saved_log_distance_ratios[subset0[0]]
         live_likelihood_index = saved_likelihood_indices[subset0[0]]
         # set the scale based on the lowest point
 
@@ -776,6 +778,7 @@ def _configure_batch_sampler(main_sampler,
         batch_sampler.scale = live_scale
         batch_sampler.live_blobs = live_blobs
         batch_sampler.distance_insertion_index = live_distance_index
+        batch_sampler.log_distance_ratio = live_log_distance_ratio
         batch_sampler.likelihood_insertion_index = live_likelihood_index
 
         batch_sampler.update_bound_if_needed(logl_min)
@@ -1111,7 +1114,7 @@ def results(self):
                 ('bound_iter', np.array(self.saved_run['bounditer'])))
             results.append(
                 ('samples_bound', np.array(self.saved_run['boundidx'])))
-            for key in ['scale', 'distance_insertion_index', 'likelihood_insertion_index']:
+            for key in ['scale', 'distance_insertion_index', 'log_distance_ratio', 'likelihood_insertion_index']:
                 results.append((key, np.array(self.saved_run[key])))
 
         return Results(results)
@@ -1367,6 +1370,7 @@ def sample_initial(self,
                             bounditer=results.bounditer,
                             scale=self.sampler.scale,
                             distance_insertion_index=self.sampler.distance_insertion_index,
+                            log_distance_ratio=self.sampler.log_distance_ratio,
                             likelihood_insertion_index=self.sampler.likelihood_insertion_index,
                             )
 
@@ -1415,6 +1419,7 @@ def sample_initial(self,
                             bounditer=results.bounditer,
                             scale=self.sampler.scale,
                             distance_insertion_index=-1,
+                            log_distance_ratio=-1,
                             likelihood_insertion_index=-1,
                             )
 
@@ -1626,6 +1631,7 @@ def sample_batch(self,
                      bounditer=results.bounditer,
                      scale=batch_sampler.scale,
                      distance_insertion_index=batch_sampler.distance_insertion_index,
+                     log_distance_ratio=batch_sampler.log_distance_ratio,
                      likelihood_insertion_index=batch_sampler.likelihood_insertion_index,
                      )
             self.new_run.append(D)
@@ -1678,6 +1684,7 @@ def sample_batch(self,
                      bounditer=results.bounditer,
                      scale=batch_sampler.scale,
                      distance_insertion_index=-1,
+                     log_distance_ratio=-1,
                      likelihood_insertion_index=-1,
                      )
             self.new_run.append(D)
@@ -1713,7 +1720,9 @@ def combine_runs(self):
         for k in [
                 'id', 'u', 'v', 'logl', 'nc', 'boundidx', 'it', 'bounditer',
                 'n', 'scale', 'blob', 'logvol',
-                'distance_insertion_index', 'likelihood_insertion_index',
+                'distance_insertion_index',
+                'log_distance_ratio',
+                'likelihood_insertion_index',
         ]:
             saved_d[k] = np.array(self.saved_run[k])
             new_d[k] = np.array(self.new_run[k])
@@ -1766,7 +1775,9 @@ def combine_runs(self):
             for k in [
                     'id', 'u', 'v', 'logl', 'nc', 'boundidx', 'it',
                     'bounditer', 'scale', 'blob',
-                    'distance_insertion_index', 'likelihood_insertion_index',
+                    'distance_insertion_index',
+                    'log_distance_ratio',
+                    'likelihood_insertion_index',
             ]:
                 add_info[k] = add_source[k][add_idx]
             self.saved_run.append(add_info)

py/dynesty/sampler.py

@@ -105,6 +105,7 @@ def __init__(self,
         # set to none just for qa
         self.scale = None
         self.distance_insertion_index = -1
+        self.log_distance_ratio = -1
         self.likelihood_insertion_index = -1
         self.method = None
         self.kwargs = {}
@@ -263,7 +264,7 @@ def results(self):
                                         dtype=int)))
             results.append(('samples_bound',
                             np.array(self.saved_run['boundidx'], dtype=int)))
-        for key in ['scale', 'distance_insertion_index', 'likelihood_insertion_index']:
+        for key in ['scale', 'distance_insertion_index', 'log_distance_ratio', 'likelihood_insertion_index']:
             results.append((key, np.array(self.saved_run.D[key])))
 
         return Results(results)
@@ -425,7 +426,18 @@ def _distance_insertion_index(self, start, point):
         norms = np.std(self.live_u, axis=0)
         distance = np.linalg.norm((point - start) / norms)
         all_distances = np.array([np.linalg.norm((start - u) / norms) for u in self.live_u])
-        return sum(all_distances < distance)
+        idx = sum(all_distances < distance)
+        if distance == 0:
+            log_ratio = np.inf
+        else:
+            other = self.live_u[np.random.choice(len(self.live_u))]
+            other_distance = np.linalg.norm((other - start) / norms)
+            while other_distance == 0:
+                other = self.live_u[np.random.choice(len(self.live_u))]
+                other_distance = np.linalg.norm((other - start) / norms)
+            alt_distance = np.linalg.norm((point - other) / norms)
+            log_ratio = self.ndim * (np.log(other_distance / distance) + np.log(other_distance / alt_distance)) / 2
+        return log_ratio, idx
 
     def _likelihood_insertion_index(self, logl):
         """
@@ -469,7 +481,7 @@ def _new_point(self, loglstar):
                 # If it's not empty we are just accumulating the
                 # the history of evaluations
                 self.update_proposal(blob, update=self.nqueue <= 0)
-                self.distance_insertion_index = self._distance_insertion_index(blob["start"], u)
+                self.log_distance_ratio, self.distance_insertion_index = self._distance_insertion_index(blob["start"], u)
             self.likelihood_insertion_index = self._likelihood_insertion_index(logl)
 
             # the reason I'm not using self.ncall is that it's updated at
@@ -602,6 +614,7 @@ def add_live_points(self):
                         scale=self.scale,
                         blob=old_blob,
                         distance_insertion_index=-1,
+                        log_distance_ratio=-1,
                         likelihood_insertion_index=-1,
                         ))
             self.eff = 100. * (self.it + i) / self.ncall  # efficiency
@@ -634,7 +647,10 @@ def _remove_live_points(self):
                 for k in [
                         'id', 'u', 'v', 'logl', 'logvol', 'logwt', 'logz',
                         'logzvar', 'h', 'nc', 'boundidx', 'it', 'bounditer',
-                        'scale', 'blob', 'distance_insertion_index', 'likelihood_insertion_index',
+                        'scale', 'blob',
+                        'distance_insertion_index',
+                        'log_distance_ratio',
+                        'likelihood_insertion_index',
                 ]:
                     del self.saved_run[k][-self.nlive:]
         else:
@@ -926,6 +942,7 @@ def sample(self,
                          scale=self.scale,
                          blob=old_blob,
                          distance_insertion_index=self.distance_insertion_index,
+                         log_distance_ratio=self.log_distance_ratio,
                          likelihood_insertion_index=self.likelihood_insertion_index,
                          ))
 

py/dynesty/utils.py

@@ -278,6 +278,7 @@ def __init__(self, dynamic=False):
             'bounditer',  # active bound at a specific iteration
             'scale',  # scale factor at each iteration
             'distance_insertion_index',
+            'log_distance_ratio',
             # number of points less distant from the starting point than the last inserted point
             'likelihood_insertion_index',
             # number of points with lower likelihood than the last inserted point
@@ -622,6 +623,11 @@ def print_fn_fallback(results,
     ('distance_insertion_index', 'array[int]',
      "The number of live points closer to the start point than "
      "the new point", 'niter'),
+    ('log_distance_ratio', 'array[float]',
+     "The log of the ratio distance of the new point to the start point to "
+     "the distance of the start point to a random other live point. This is "
+     "used as a diagnostic for independence of the new sample",
+     'niter'),
     ('likelihood_insertion_index', 'array[int]',
      "The number of live points with likelihood less than "
      "the new point", 'niter'),