Relbin resolution diagnostic - #5409
Open
ahnitz wants to merge 13 commits into
Open
Conversation
Relative binning assumes the ratio of the waveform to the fiducial one is linear across a frequency bin, but nothing checked that the bins chosen by epsilon actually deliver that, so a value too coarse for the problem passed silently. Evaluate the ratio at the bin midpoints for a few draws from the prior and compare it to the linear model the approximation assumes. This measures the approximation directly, and costs two sparse waveforms per draw: on GW170817 the check takes the time of about 29 likelihood evaluations, against the millions an analysis uses. Replaces max_curvature_from_reference, which was never called. That estimated the same thing from a second difference of the ratio over the bin edges, but the edges are not evenly spaced, so the estimate kept a term in the first derivative and did not fall reliably as bins were added.
The resolution diagnostic generates a waveform of its own at the bin midpoints, and passed the model's current parameters straight to the generator. Those need not be single numbers: a marginalized parameter is held as the whole vector of samples drawn for it, and some samplers hand a batch of points in at once. The generator takes one point at a time and raised a TypeError on the vector rather than saying so, which took out epsilon='auto' wherever such a parameter reached the check. Take the first of each. The error measured is relative to the size of the ratio it is an error in, so a parameter that only scales the waveform, distance above all, divides back out, and the sky and time parameters never reach the generator.
check_bin_resolution runs at the end of every Relative construction and draws a few points from the prior. A model built without one is given a stand-in object rather than None, so the guard against there being no prior did not catch it and construction raised. Ask whether the prior can be drawn from rather than whether it is there.
Of no use to someone reading the test for the first time; the history is in the commit that fixed it.
Three corrections to the resolution check, none of which change the number it reports. The linear model the check compares against was written out by hand as a weight and two multiply-adds. numpy.interp does exactly that and takes complex values for the ordinates, so it says what the line means and is faster everywhere that matters: at the 62 bins this model actually has, 0.5 us against 3.2 us, and 248 us against 639 us at 1e5 bins. The two agree to 3.6e-18 absolute, 1.6e-16 relative, on both detectors of the test fixture, which is rounding. The docstring said the check costs two waveforms per draw. Counting the calls into the generator, it is three for the two-detector fixture: one at the bin edges, from the likelihood's own call, which is shared across detectors whose edges coincide, plus one at the midpoints for each detector. Ten draws measure at 25 times the wall time of a likelihood evaluation, which is what the earlier "about 29 likelihood evaluations" was reporting and is consistent with three per draw, not two. The GW170817 numbers quoted for the threshold imply factors of 32, 24, 25 and 57, not "roughly 25". Quote the four and say the factor is not a constant, which the next sentence about signal to noise already implied. test_check_is_cheap asserted only that the returned value was finite and non-negative, so it measured none of the cost its docstring claimed. It now counts the calls into the waveform generator and the number of sample points asked of each. Generating at the data's frequencies rather than the bin midpoints fails it, as does one extra waveform per draw; both were checked by mutation. Also itertools.pairwise for the two adjacent-pair loops, and the module docstring said the file runs on cached data when it builds its own simulated noise. It runs in 1.8 seconds.
On this branch interpolation_error_from_reference has one caller, check_bin_resolution, which evaluates a single prior draw at a time. Every parameter it passes is one number, so reducing them was dead code. The caller that can pass a vector is the signal-to-noise weighted screen, which runs while the sampler does and belongs to the branch that adds it. The reduction goes there, with the test for it.
The diagnostic is interpolation_error_from_reference: it evaluates the waveform ratio at the bin midpoints and returns how far it departs from the linear model the bins assume. What it does not do any more is run itself at construction over a set of prior draws. A prior of any size draws mostly from where the posterior is not, so the error it found there was not the error anywhere the sampler would go; and the misconfiguration it warned about does not go unnoticed anyway. With the fiducial mass moved from 1.4 to 1.39, the log likelihood ratio at the true parameters falls from 52.6 to -12.9, and at 1.37 to -291. Nobody misses that, and a line in a log adds nothing to it. What the number is good for is comparing one binning against another, which is what the caller can now do with it.
Relative binning assumes the ratio to the fiducial waveform is linear between bin edges. Replace max_curvature_from_reference, which reported a second difference in units that do not correspond to anything, with a direct measure of that assumption: predict each interior edge from its two neighbours and compare with the value the likelihood already evaluated. The result is relative to the size of the ratio, so it is comparable across detectors and bins. Reaching across two bins, it reports the error of a binning twice as coarse, about four times the error within a bin. No waveform is generated to do this, so it costs nothing beyond the likelihood call and works on both the polarization and the detector response waveform paths.
The interpolation error depends on where in the parameter space a call lands, so a single measurement does not stand in for a run. Add check_interpolation_error, off by default, which measures it on every likelihood call and keeps the largest in max_interpolation_error, written to the output file. It is a debugging aid for choosing epsilon. It costs a fixed 0.02-0.03 ms, independent of the marginalization: 22% of an unmarginalized likelihood call, 0.3% of one marginalized over time and sky with 1000 samples. That is why it is opt-in. Setting self.wf_ret was missing on the detector response path, which had gone unnoticed because its only reader was unused. Give get_waveforms one exit so both paths record it.
check_interpolation_error described the mechanism rather than the thing being checked. The question the option answers is whether the bins are fine enough for the heterodyne, which is the word relbin.py and the model documentation already use. check_interpolation_error -> check_heterodyne_bins max_interpolation_error -> max_heterodyne_error interpolation_error_from_reference -> heterodyne_bin_error
The model is copied into each worker, so a maximum kept on it covers only whatever calls that copy happened to make, and what reaches the output file is whichever copy wrote last. Log each measurement instead and leave the reader to take the maximum over the run. This drops max_heterodyne_error and the attribute written to the file.
They needed a fixture that is already there. TestModels sets up the reference GW170817 analysis with the same static parameters, priors and fiducial waveform these tests were building for themselves, so reusing it drops the whole fixture and tests against real data rather than simulated noise. This removes test_heterodyne_bins.py.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This checks (optionally) the accuracy of the heterodyne frequency bins by comparing the predicted ratios from those at a 2x courser resolution. This is done from the info on hand when a likelihood is calculation. For now it is just a piece of debugging information that can be sent to logging with an option to the model.
This is prestaging some work later though to hopefully use this to set the binning dynamically.