Repository navigation
Convolve with a filter bank in peak finding - #70
Merged
Merged
Conversation
Applying N psfs to one image with N filter structures computes the forward FFT of that image N times, and all N results are identical. A bank keeps one pair of FFTW plans and N precomputed psf transforms, transforms the image once, and then does N multiply-and-invert passes. The transform count per call goes from 2N to N+1. image_fft is separate from work_fft because convolve() does its complex multiply in place. Reusing one buffer would destroy the very transform the bank exists to share. Measured on 270x270, the size the 256x256 test data pads to, interleaving the two arms so that neither gets a whole block of warm cache: N separate bank speedup 2N/(N+1) 1 0.568 ms 0.574 ms 0.99x 1.00x 2 1.267 0.935 1.36x 1.33x 3 1.911 1.317 1.45x 1.50x 5 3.339 1.997 1.67x 1.67x 8 5.499 3.106 1.77x 1.78x So it tracks the predicted ratio, and at N = 1 it correctly does nothing. There is no point using a bank for a single psf. The bank does not memoize. convolveMemo() only helps when the same filter sees nearly the same image twice in a row, and instrumenting the spliner peak finder showed that never happening once in 220 calls - the image changes every iteration as peaks are subtracted. Nothing is given up by leaving it out. test_matched_filter_bank1 asserts array_equal against N separate filters rather than allclose, since the bank changes only where the forward transform happens, not the arithmetic. Also lists the existing tests in the __main__ block, which called only test_matched_filter5.
PeakFinderArbitraryPSF convolves at each of the z values in z_value, and both convolution inputs are the same at every one of them: bg_var for the variance filters, and image - background - fit_peaks for the foreground filters. Each z value was transforming those same two arrays again. With one bank per filter type the two inputs are transformed once per call instead of once per z value, so a run with N z values goes from 4N transforms per peak finding pass to 2N + 2. The foreground subtraction was also being recomputed inside the loop to produce the same array every time, so it is now done once where the convolution happens. Measured on test_spliner.dax with test_spliner_dh.xml, which uses five z values, timing peakFinder() only and alternating the two versions so that neither arm gets a whole block of warm cache: bank 0.176s master 0.235s 1.33x bank 0.170s master 0.237s 1.39x bank 0.172s master 0.291s 1.69x bank 0.219s master 0.291s 1.33x The convolutions themselves go 1.67x at N = 5; peak finding as a whole gains less because findMaxima and the masking are unchanged. That works out to roughly a sixth of steady state analysis time for this configuration. At the default z_value of [0.0] there is one filter per bank and nothing changes, which is the expected result rather than a disappointment - the redundancy only exists when there is more than one z value. Verification. Analyzing with fftw_estimate = 1 gives a byte for byte identical HDF5 before and after, across all 234 datasets. That qualifier matters. With the default fftw_estimate = 0 the output is NOT reproducible, on master or on this branch: two master runs of the same file differ in 2 to 8 datasets at up to 1.1e-15 relative, which is larger than the 5.0e-16 that separates master from this branch. FFTW_MEASURE picks its algorithm from plan time benchmarks, and that choice varies with timing noise. Pinning the planner removes the variation from both arms and leaves them identical, which is the real statement that this change is arithmetically inert. Full suite green at 280 tests.
Same redundancy as the spliner peak finder, one dimension over. The filters were stored indexed by z value and then by channel, and the loops ran z on the outside, so each channel's image was transformed once per z value even though the image does not depend on z. The filters are now one bank per channel, each holding that channel's psf at every z value, and each channel's image is convolved once per pass. For N z values and C channels that takes the transform count from 4NC per pass to 2C(N + 1). The loops still accumulate z on the outside and channel on the inside, so the order the per-channel contributions are summed in is unchanged. That matters because those sums are floating point and reordering them would perturb the result. MPPeakFinderDao is left on the same code path but has a single z value, so its banks hold one psf each and it does exactly what it did before. That is expected: like 3D-DAOSTORM, it convolves each image once, and there is no shared transform to hoist. Verified against the multiplane_spline diagnostic, two channels and three z values, with fftw_estimate = 1 so that the FFTW planner cannot vary between runs: all 191 datasets in the output HDF5 are identical before and after. Timing the peak finding step only, alternating the arms: bank 0.539s master 0.699s 1.30x bank 0.635s master 0.684s 1.08x bank 0.627s master 0.676s 1.08x Less than the 1.5x the convolutions themselves get at three z values, because multiplane peak finding also does the affine transforms into the channel 0 frame and the accumulation, neither of which changed. Full suite green at 282 tests.
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
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.
Peak finding convolves the image with the PSF at each of the z values in
z_value. The image does not depend on z, so every one of those convolutionscomputed the same forward FFT of the same array and threw it away. A filter
bank keeps one pair of FFTW plans and N precomputed PSF transforms, transforms
the image once, and then does N multiply-and-invert passes.
For N z values that takes the peak finding pass from 4N transforms to 2N + 2.
Where this helps, and where it does not
It only helps where several filters are applied to the same image in one
pass. That is the spliner z loop and the multiplane z/channel loop.
3D-DAOSTORM gets nothing from this, by construction. Its three filters
(background, foreground, foreground variance) each act on a different image, so
there is no shared transform to hoist. The same is true of
MPPeakFinderDao,which has a single z value, and of spliner at the default
z_valueof[0.0].Those paths run the new code with single-PSF banks and do exactly what they did
before. That is the expected result rather than a disappointment - the
redundancy only exists when there is more than one z value.
Measurements
The bank itself, at 270x270, which is what 256x256 data pads to. Arms
interleaved so neither gets a block of warm cache:
End to end, timing the peak finding step only and alternating the two versions:
test_spliner_dh.xml, five z values: 1.33x, 1.39x, 1.69x, 1.33x.Roughly a sixth of steady state analysis time.
1.08x. Less than the 1.5x the convolutions get, because multiplane peak
finding also does the affine transforms into the channel 0 frame and the
accumulation, and neither changed.
Verification
Both call sites produce byte for byte identical output:
multiplane_spline diagnostic since there is no end to end multiplane test.
Both comparisons require
fftw_estimate = 1, and that qualifier is worthstating plainly, because it turned up something about the existing code:
with the default
fftw_estimate = 0this project's output is notreproducible run to run, on master or on this branch. Two master runs of the
same file differ in 2 to 8 datasets at up to 1.1e-15 relative, which is larger
than the 5.0e-16 that separated master from this branch before the planner was
pinned.
FFTW_MEASUREchooses its algorithm from plan time benchmarks and thatchoice moves with timing noise. Pinning the planner removes the variation from
both arms and leaves them identical, which is the real statement that this
change is arithmetically inert.
The practical consequence is that "re-run and diff" is not a valid regression
check for anything downstream of these filters unless the planner is pinned
first.
To be clear about the default itself:
FFTW_MEASUREearns its planning costseveral times over and should stay. Measured per convolution, in separate
processes so that the two arms cannot share FFTW wisdom, at the padded sizes
fftFriendlyMargin()actually produces:On a 10k frame movie that is half a minute to six minutes saved for under a
second of planning.
One behaviour change
The banks do not memoize.
convolveMemo()only helps when a filter sees nearlythe same image twice in a row, and instrumenting the spliner peak finder showed
that never happening once in 220 calls, because the image changes every
iteration as peaks are subtracted. Nothing measurable is given up, but it is a
removal rather than a pure refactor.
Full suite green at 282 tests, the two new ones covering the bank against
separate filters with
array_equalrather thanallclose.