Let relative binning choose its own frequency resolution

epsilon had to be picked by hand with nothing to pick it from: the
default in code is 0.5 while the examples use 0.03 to 0.1, a factor of
sixteen with no guidance on which suits a given analysis.

Accept epsilon = auto, which halves epsilon until the bin resolution
check passes, so the resolution is set by the problem. On GW170817 that
settles on 0.5 and stops; with a fiducial waveform far from the signal it
keeps refining and then says so, since bins cannot rescue a bad fiducial.

The bin layout moves into its own method so it can be redone, which also
lets the shifted data used to build it be dropped once it has been. The
epsilon that was used is now written to the output file, which the
metadata docstring already claimed.

Author: ahnitz

pycbc/inference/models/relbin.py

@@ -158,9 +158,11 @@ class Relative(DistMarg, BaseGaussianNoise):
         an intrinsic parameter compatible with the chosen approximant.
     gammas : array of floats, optional
         Frequency powerlaw indices to be used in computing frequency bins.
-    epsilon : float, optional
+    epsilon : float or 'auto', optional
         Tuning parameter used in calculating the frequency bins. Lower values
-        will result in higher resolution and more bins.
+        will result in higher resolution and more bins. If 'auto', bins are
+        added until they resolve the waveform ratio, so that the resolution
+        is set by the problem rather than by hand.
     earth_rotation: boolean, optional
         Default is False. If True, then vary the fp/fc polarization values
         as a function of frequency bin, using a predetermined PN approximation
@@ -220,6 +222,10 @@ def __init__(
             if self.fid_params[k] == 'REPLACE':
                self.fid_params.pop(k)
 
+        # the shifted data and frequency limits the bin layout is built
+        # from, kept only until the layout has been settled on
+        shifted, limits = {}, {}
+
         for ifo in data:
             # store data and frequencies
             d0 = self.data[ifo]
@@ -291,22 +297,83 @@ def __init__(
             self.h00[ifo] = numpy.array(curr_wav) # * tshift
             data_shifted = self.data[ifo] * numpy.conjugate(tshift)
 
+            shifted[ifo] = data_shifted
+            limits[ifo] = (f_lo, f_hi)
+
+        layout = (shifted, limits, gammas, earth_rotation,
+                  int(earth_rotation_mode))
+        if str(epsilon).lower() == 'auto':
+            self.epsilon = self.refine_epsilon(*layout)
+        else:
+            self.epsilon = float(epsilon)
+            self.setup_bin_layout(self.epsilon, *layout)
+            self.check_bin_resolution()
+
+    def setup_bin_layout(self, epsilon, shifted, limits, gammas,
+                         earth_rotation, earth_rotation_mode):
+        """Place the frequency bins and compute the summary data for them.
+        """
+        for ifo in self.data:
             logging.info("Computing frequency bins")
+            f_lo, f_hi = limits[ifo]
             fbin_ind = setup_bins(
                 f_full=self.f[ifo], f_lo=f_lo, f_hi=f_hi,
-                gammas=gammas, eps=float(epsilon),
+                gammas=gammas, eps=epsilon,
             )
             logging.info("Using %s bins for this model", len(fbin_ind))
 
             self.fedges[ifo] = self.f[ifo][fbin_ind]
             self.edges[ifo] = fbin_ind
-            self.init_from_frequencies(data_shifted, self.h00, fbin_ind, ifo)
+            self.init_from_frequencies(shifted[ifo], self.h00, fbin_ind, ifo)
             self.antenna_time[ifo] = self.setup_antenna(
                                         earth_rotation,
-                                        int(earth_rotation_mode),
+                                        earth_rotation_mode,
                                         self.fedges[ifo])
         self.combine_layout()
-        self.check_bin_resolution()
+
+    def refine_epsilon(self, *layout, start=1.0, threshold=1e-3, smallest=0.01):
+        """Halve epsilon until the bins resolve the waveform ratio.
+
+        Bins are added only where the model cannot do without them, so a
+        problem needing few bins does not pay for one that needs many.
+        Each step costs a bin layout and a handful of sparse waveforms,
+        which is small against the analysis that follows.
+
+        Parameters
+        ----------
+        start : float, optional
+            The first, coarsest value to try.
+        threshold : float, optional
+            Stop once the interpolation error falls below this. See
+            :py:meth:`check_bin_resolution` for what it means.
+        smallest : float, optional
+            Give up at this value rather than refining without end.
+
+        Returns
+        -------
+        float
+            The epsilon that was settled on.
+        """
+        epsilon = start
+        while True:
+            self.setup_bin_layout(epsilon, *layout)
+            # inf so that only the final value is reported to the user
+            error = self.check_bin_resolution(threshold=numpy.inf)
+            if error <= threshold or epsilon <= smallest:
+                break
+            epsilon /= 2.
+
+        if error > threshold:
+            logging.warning(
+                "Refining epsilon stopped at %.3g without resolving the "
+                "waveform ratio, which is still off by %.3g against a "
+                "threshold of %.3g. The fiducial waveform is probably far "
+                "from the posterior.", epsilon, error, threshold)
+        else:
+            logging.info("Chose epsilon %.3g, giving %s bins and an "
+                         "interpolation error of %.3g", epsilon,
+                         len(self.fedges[list(self.data)[0]]), error)
+        return epsilon
 
     def init_from_frequencies(self, data, h00, fbin_ind, ifo):
         bins = numpy.array(
@@ -636,6 +703,7 @@ def write_metadata(self, fp, group=None):
             attrs = fp[group].attrs
         for p, v in self.fid_params.items():
             attrs["{}_ref".format(p)] = v
+        attrs["epsilon"] = self.epsilon
 
     def interpolation_error_from_reference(self):
         """ Return the largest error made by interpolating the waveform

test/test_relbin_resolution.py

@@ -127,6 +127,39 @@ def test_check_is_cheap(self):
         self.assertTrue(numpy.isfinite(value))
         self.assertGreaterEqual(value, 0.)
 
+    # 'auto' picks a resolution that suits the problem
+
+    def test_auto_meets_the_threshold(self):
+        model = self.model('auto')
+        self.assertLess(model.check_bin_resolution(), 1e-3)
+
+    def test_auto_agrees_with_a_finer_setting(self):
+        """Stopping as soon as the check passes must not cost accuracy."""
+        fine = self.model(0.005)
+        fine.update(**self.q)
+
+        auto = self.model('auto')
+        auto.update(**self.q)
+        self.assertAlmostEqual(auto.loglr, fine.loglr, delta=0.05)
+
+    def test_auto_spends_bins_only_where_needed(self):
+        """A fiducial further from the signal should buy more bins."""
+        near = self.model('auto')
+        offset = dict(self.static)
+        offset['mass1'], offset['mass2'] = 1.39, 1.36
+        far = self.model('auto', static=offset)
+
+        self.assertGreater(far.epsilon, 0.)
+        self.assertLess(far.epsilon, near.epsilon)
+        self.assertGreater(len(far.fedges['H1']), len(near.fedges['H1']))
+
+    def test_epsilon_is_recorded(self):
+        """The value used has to be recoverable, since it is now chosen."""
+        self.assertEqual(self.model(0.25).epsilon, 0.25)
+        # config files hand every option over as a string
+        self.assertEqual(self.model('0.25').epsilon, 0.25)
+        self.assertIsInstance(self.model('auto').epsilon, float)
+
 
 suite = unittest.TestSuite()
 suite.addTest(