Allocate the rate and the draws together, not the rate first

Asking for both 'auto' was the most expensive thing a user could ask
for: on a 2.2+1.6 injection at the default accuracy it cost 29.6 ms a
likelihood call against 6.6 ms for a hand-picked 8192 Hz with the draws
chosen, for the same accuracy. That is backwards, and both 'auto' is what
a user will write.

The cause was the order. choose_sample_rate ran first and spent the whole
budget on the rate, since its target is 1/(vsamples*accuracy**2) with
vsamples still at its starting value; by the time the draws were sized
the budget was already met and the cheap lever had nothing left to buy.
The rate is the expensive one - the signal to noise series is rebuilt
every call, so four times the rate costs about four times as much, while
four times the draws cost tens of a percent - so it should be the lever
of last resort, not of first.

The accuracy depends only on the product of the two, which is what makes
them exchangeable. So when the number of draws is being chosen, the rate
is asked only for a resolution the product law can be trusted at and the
draws are sized to the rest of the budget. Only if the budget needs more
draws than there is room for - MOST_VSAMPLES, or the precalculated pool,
which is the real ceiling - does the rate climb, to the resolution the
most draws available would need, capped at FLOOR_RESOLVED where the rate
stops buying anything at all. That target is either reached or it is not
and the next one would be the same number, so the iteration goes round at
most twice and then says plainly that the accuracy is out of reach and
that a bigger pool, not a higher rate, is the way out.

The resolution to stop at is measured rather than picked. Sizing the
draws by the law and measuring what comes out: at one sample across the
peak the scatter is 2.8 times what was asked, at two samples - the hard
floor, where the answer stops being biased - it is still 1.5 times, and
at four samples it is 0.99 and 0.82 at two budgets a factor of four
apart. Four is where the law starts delivering what it promises, and on
this signal it is also where the total cost of a call bottoms out, since
the draws needed rise as the rate falls. Hence SANE_RESOLVED = 4.

Both 'auto' at the default accuracy now picks 8192 Hz and 13942 draws:
29.59 -> 6.14 ms a call, 4.8x cheaper, with the scatter improving from
0.0058 to 0.0053 against the 0.005 asked for. Where the draws are given
explicitly nothing changes: the rate still buys the whole budget on its
own, by the same code and with the same message.

Author: ahnitz

pycbc/inference/models/relbin.py

@@ -69,6 +69,18 @@
 FEWEST_VSAMPLES = 32
 MOST_VSAMPLES = 32768
 
+# The least resolution to ask the sample rate for when the number of draws
+# is being chosen too. The accuracy depends only on the product of the two,
+# so the rate is wanted no finer than the point where that product law can
+# be relied on. Measured, it cannot be relied on at the hard floor of two
+# samples across the peak: with the draws sized by the law, two samples
+# deliver a scatter half again the size asked for and one sample nearly
+# three times it, while four samples deliver what they promise, at two
+# accuracy budgets a factor of four apart. Four is also where the total
+# cost of a call bottoms out, the draws needed rising as the rate falls, so
+# nothing is given up by stopping there.
+SANE_RESOLVED = 4.0
+
 
 def setup_bins(f_full, f_lo, f_hi, chi=1.0,
                eps=0.1, gammas=None,
@@ -840,6 +852,19 @@ def choose_sample_rate(self, accuracy=None, most=262144.0, **kwargs):
         neutron star, and above a few kHz the series is already most of
         what the model does.
 
+        Where the number of marginalization draws is being chosen as well,
+        the two are allocated together, because the accuracy depends only
+        on their product and they are not the same price. Four times the
+        rate costs about four times as much, while four times the draws
+        cost tens of a percent, so the cheapest way to a given product is
+        as much of it as possible from the draws. The rate is then asked
+        for no more than ``SANE_RESOLVED``, the resolution below which the
+        product law stops delivering what it promises, and the draws are
+        sized to the rest of the budget. Only if there is not room for
+        enough draws does the rate climb further, which it does until
+        either the budget is met or the resolution stops being what limits
+        the answer.
+
         Parameters
         ----------
         accuracy : float, optional
@@ -857,24 +882,74 @@ def choose_sample_rate(self, accuracy=None, most=262144.0, **kwargs):
         if accuracy is None:
             accuracy = self.marginalization_accuracy
 
-        # invert marginalization_error for the resolution wanted. Two
-        # samples is a floor rather than an accuracy: below it the peak
-        # falls between grid points and the answer is biased, not merely
-        # noisy, so no scatter budget makes that acceptable.
-        target = max(2.0, 1.0 / (self.vsamples * accuracy ** 2))
-
-        resolved = self.resolved_samples(self.ref_snr)
-        while resolved < target and self.sample_rate < most:
-            if resolved >= 2.0:
-                # resolved well enough to say how much more is needed
-                steps = max(1, int(numpy.ceil(numpy.log2(target / resolved))))
-            else:
-                steps = 1
-            self.set_sample_rate(min(self.sample_rate * 2 ** steps, most),
-                                 **kwargs)
+        if self.adapt_vsamples:
+            # the draws will buy the accuracy, so ask the rate only for a
+            # resolution the product law can be trusted at
+            ceiling = self.samples_ceiling(**kwargs)
+            target = SANE_RESOLVED
+        else:
+            # invert marginalization_error for the resolution wanted. Two
+            # samples is a floor rather than an accuracy: below it the peak
+            # falls between grid points and the answer is biased, not merely
+            # noisy, so no scatter budget makes that acceptable.
+            ceiling = None
+            target = max(2.0, 1.0 / (self.vsamples * accuracy ** 2))
+
+        while True:
             resolved = self.resolved_samples(self.ref_snr)
+            while resolved < target and self.sample_rate < most:
+                if resolved >= 2.0:
+                    # resolved well enough to say how much more is needed
+                    steps = max(1,
+                                int(numpy.ceil(numpy.log2(target / resolved))))
+                else:
+                    steps = 1
+                self.set_sample_rate(min(self.sample_rate * 2 ** steps, most),
+                                     **kwargs)
+                resolved = self.resolved_samples(self.ref_snr)
+
+            if ceiling is None:
+                break
+
+            # what the draws would have to be to meet the budget here. If
+            # there is no room for that many the rate has to make up the
+            # difference after all: ask it for the resolution the most draws
+            # available would need, and no more than the resolution past
+            # which the rate buys nothing whatever it is asked. The new
+            # target is reached or it is not, and either way the one after
+            # it would be the same number, so this goes round at most twice.
+            wanted = 1.0 / (accuracy ** 2 * min(resolved, FLOOR_RESOLVED))
+            if wanted <= ceiling:
+                break
+            harder = min(1.0 / (accuracy ** 2 * ceiling), FLOOR_RESOLVED)
+            if harder <= target:
+                break
+            target = harder
 
         rate = self.sample_rate
+        if ceiling is not None:
+            if wanted > ceiling:
+                logging.warning("Could not reach the accuracy asked of the "
+                                "marginalization: %.1f samples across the "
+                                "peak at %s Hz would need %.3g draws for a "
+                                "scatter of %.3g and there is room for %s, "
+                                "which leaves about %.3g. A larger "
+                                "precalculate_marginalization_points, or a "
+                                "looser marginalization_accuracy, is the way "
+                                "out; a higher rate is not",
+                                resolved, rate, wanted, accuracy, ceiling,
+                                1.0 / (ceiling
+                                       * min(resolved,
+                                             FLOOR_RESOLVED)) ** 0.5)
+            else:
+                logging.info("Chose a sample rate of %s Hz, spreading the "
+                             "peak over %.1f samples: enough resolution for "
+                             "the draws to be trusted, and the %.3g asked "
+                             "for is bought with the draws from there, which "
+                             "is much the cheaper of the two",
+                             rate, resolved, accuracy)
+            return rate
+
         error = self.marginalization_error(resolved)
         if resolved < target:
             logging.warning("Could not resolve the peak of the likelihood in "
@@ -906,6 +981,20 @@ def choose_sample_rate(self, accuracy=None, most=262144.0, **kwargs):
                          rate, resolved, error, accuracy)
         return rate
 
+    def samples_ceiling(self, precalculate_marginalization_points=False,
+                        **kwargs):
+        """The most marginalization draws that could be asked for.
+
+        Taken from the model's keywords rather than from what has been built
+        so far, because the rate is chosen before the pool of points exists
+        and how many draws there is room for is what decides whether the
+        rate has to climb.
+        """
+        if precalculate_marginalization_points:
+            return min(MOST_VSAMPLES,
+                       int(float(precalculate_marginalization_points)))
+        return MOST_VSAMPLES
+
     def wanted_ess(self, accuracy=None):
         """How many effective marginalization samples the accuracy asks for.
 

test/test_marg_vsamples.py

@@ -33,6 +33,7 @@
 """
 
 import copy
+import time
 import unittest
 
 import numpy
@@ -307,36 +308,106 @@ def test_points_far_from_the_best_do_not_move_it(self):
         self.assertLess(float(model.loglr), model.marg_best_loglr - 20.)
         self.assertEqual(model.vsamples, settled)
 
-    def test_both_knobs_on_auto_reach_the_accuracy_together(self):
-        """The configuration someone would actually write.
-
-        The resolution is chosen once at setup, from the number of samples
-        it starts with; the number of samples is then chosen for the
-        resolution that was settled on, and followed from there. The two
-        buy the same accuracy at very different prices, so what matters is
-        that between them they deliver it.
-        """
-        accuracy = 0.03
+    def auto_rate_model(self, vsamples='auto', accuracy=0.01):
+        """The same model with the rate left to the model as well."""
         variable = ['distance', 'inclination', 'tc']
         prior = JointDistribution(
             list(variable), SinAngle(inclination=None),
             Uniform(distance=(10, 300)),
             Uniform(tc=(TC - HALFWIDTH, TC + HALFWIDTH)))
-        model = models.RelativeTimeDom(
+        return models.RelativeTimeDom(
             list(variable), copy.deepcopy(self.data),
             low_frequency_cutoff=self.flow, psds=self.psds,
             static_params=self.static, prior=prior,
             fiducial_params={'mass1': INJ['mass1'], 'mass2': INJ['mass2'],
                              'tc': TC},
             epsilon=0.1, marginalize_vector_params='tc', sample_rate='auto',
             marginalization_accuracy=accuracy,
-            marginalize_vector_samples='auto')
+            marginalize_vector_samples=vsamples)
+
+    def cost(self, model, ncall=40):
+        """Seconds a likelihood evaluation takes."""
+        numpy.random.seed(31)
+        start = time.time()
+        for _ in range(ncall):
+            model.update(**self.point)
+            model.loglr
+        return (time.time() - start) / ncall
+
+    def test_asking_for_both_does_not_buy_the_expensive_one(self):
+        """The accuracy depends only on the product of resolution and
+        draws, and the two are not the same price: the signal to noise
+        series is rebuilt every call, so the rate is close to linear in
+        cost, while quadrupling the draws costs tens of a percent. Asking
+        for both must therefore not cost more than choosing a low rate by
+        hand and asking only for the draws -- which is what happens if the
+        rate is allowed to spend the whole budget before the draws are
+        sized. A factor of two is allowed on a timing taken in one process
+        on one signal; the point is that it is not the fourfold or more
+        that leaning on the rate costs.
+        """
+        both = self.auto_rate_model()
+        byhand = self.model(accuracy=0.01)
+        self.drive(both, 200)
+        self.drive(byhand, 200)
+        together, apart = self.cost(both), self.cost(byhand)
+        self.assertLess(together, 2.0 * apart,
+                        "%s Hz and %s samples took %.2f ms against %.2f for "
+                        "%s Hz and %s samples"
+                        % (both.sample_rate, both.vsamples, together * 1e3,
+                           apart * 1e3, byhand.sample_rate, byhand.vsamples))
+
+    def test_asking_for_both_still_reaches_the_accuracy(self):
+        """Spending less on the rate is only worth having if the accuracy
+        still arrives, so the scatter is measured rather than predicted.
+
+        This is also the configuration someone would actually write: the
+        rate settled once at setup, the draws sized for it and then
+        followed, and the two together delivering what was asked.
+        """
+        accuracy = 0.01
+        model = self.auto_rate_model(accuracy=accuracy)
         trail = self.drive(model, 400)
         self.assertEqual(trail[-100:], [trail[-1]] * 100)
         self.assertLess(self.scatter(model), 1.5 * accuracy,
                         "%s Hz and %s samples"
                         % (model.sample_rate, model.vsamples))
 
+    def test_an_explicit_count_leaves_the_rate_choice_alone(self):
+        """A number of draws given explicitly is a number the model may not
+        change, so the rate is the only thing left to buy the accuracy with
+        and it has to buy all of it, exactly as before this option existed.
+        Asking for the draws as well must then cost strictly less rate.
+        """
+        accuracy = 0.01
+        explicit = self.auto_rate_model(vsamples=1000, accuracy=accuracy)
+        # the rate alone, against the resolution the accuracy law asks of it
+        self.assertGreaterEqual(
+            explicit.resolved_samples(explicit.ref_snr),
+            1.0 / (1000 * accuracy ** 2))
+        chosen = self.auto_rate_model(accuracy=accuracy)
+        self.assertLess(chosen.sample_rate, explicit.sample_rate)
+
+    def test_the_rate_climbs_only_when_the_draws_run_out_of_room(self):
+        """The rate is the lever of last resort.
+
+        While there is room for the draws the resolution asked for is the
+        least that can be trusted; when the budget needs more draws than
+        there can be, the rate is the only thing left and it climbs. Which
+        way round it went is visible in the draws: they are up against their
+        ceiling in the second case and nowhere near it in the first.
+        """
+        from pycbc.inference.models.relbin import (MOST_VSAMPLES,
+                                                   SANE_RESOLVED)
+        loose = self.auto_rate_model(accuracy=0.01)
+        self.assertLess(loose.resolved, 2.0 * SANE_RESOLVED)
+        self.assertLess(loose.vsamples, MOST_VSAMPLES / 2)
+
+        tight = self.auto_rate_model(accuracy=0.002)
+        self.assertGreater(tight.resolved, SANE_RESOLVED)
+        self.assertGreater(tight.sample_rate, loose.sample_rate)
+        self.assertGreater(tight.vsamples, MOST_VSAMPLES / 2)
+
     def test_an_explicit_number_is_left_alone(self):
         """Asking for a number of samples must still get exactly that."""
         model = self.model(vsamples=777)