[ESSDIFFRACTION] Align powder focussing 2theta bins with requested TwoThetaBins - #725
[ESSDIFFRACTION] Align powder focussing 2theta bins with requested TwoThetaBins#725SimonHeybrock wants to merge 7 commits into
Conversation
Focussing binned into 1024 uniform 2theta bins spanning the detector range, independently of the 2theta binning requested for the result. Grouping into the requested bins then assigned each focussing bin as a whole to the bin containing its center, so the number of bins per group varied periodically. For BEER this produced spikes of the order of the signal itself. The focussing grid is now built as a refinement of TwoThetaBins: every requested edge occurs exactly in the grid, and each requested bin is subdivided until the sub-bins reach the resolution needed to reconstruct wavelength for monitor normalization. Grouping is then an exact sum over whole sub-bins. Outside the requested range the grid continues at the base resolution so that integrating over 2theta still covers the full detector. TwoThetaBins gains a None default, meaning the data will not be grouped by 2theta, so workflows that only produce 1D d-spacing output are unaffected. Both focussing and grouping now convert between the units of the requested edges and of the two_theta coordinate. Requesting bins in degrees against a coordinate in radians previously raised. Fixes #670
github-actions
Bot
changed the title
The number only matters for histogrammed data normalized by a monitor histogram: the monitor is looked up once per (d-spacing, two-theta) bin and the lookup is piecewise constant on the monitor bins. Data that keeps its events carries its own wavelength and is unaffected.
The bins are recomputed for every chunk in a stream processor, since they depend on the two-theta coordinate of the chunked detector data. Subdividing the requested bins one at a time cost 1.6 ms for 180 requested bins, a third of the cost of the histogramming step it feeds. Computing the sub-bin positions in one go makes the cost independent of the number of requested bins, at around 0.4 ms, half of which is the min/max over the two-theta coordinate.
scipp::bin and scipp::hist have a fast path for linearly spaced bin edges. Continuing the grid beyond the requested range at the base width rather than at the width of the sub-bins broke that, costing 12% of the histogramming step. Two-theta is a per-pixel coordinate, so the binary search is amortized over pixels rather than events; the same loss is a factor of 8-10 where the binning coordinate is an event coordinate. Extending at the width of the outermost sub-bins restores linear spacing whenever the requested bins are linearly spaced. A sub-bin is never narrower than half the base width, so this at most doubles the number of bins outside the requested range.
The docstring stated that scipp can exploit linear spacing without saying what for.
jl-wynen
left a comment
jl-wynen
left a comment
There was a problem hiding this comment.
This seems to work. The main issue I see is that it creates a dependency between focussing and grouping. When users call these functions directly, they might use incompatible bins.
|
|
||
| def _extend_edges(start: float, step: float, limit: float) -> np.ndarray: | ||
| """Ascending edges reaching from ``start`` (exclusive) past ``limit``.""" | ||
| n = max(int(np.ceil((limit - start) / step)), 0) |
There was a problem hiding this comment.
What is max for? How can the argument be negative?
There was a problem hiding this comment.
It can be negative: it happens when the requested bins extend past the detector range, e.g. TwoThetaBins covering 0-180 deg on a detector spanning 10-170 deg, so the limit is already reached and no extension edges are wanted. But np.arange(1, n + 1) is empty for n <= 0 anyway, so max was doing nothing. Dropped it and named the case in the docstring instead. There is now a test with bins reaching beyond the detector, which also pins that the grid still covers the full range.
| """ | ||
|
|
||
|
|
||
| def _extend_edges(start: float, step: float, limit: float) -> np.ndarray: |
There was a problem hiding this comment.
| def _extend_edges(start: float, step: float, limit: float) -> np.ndarray: | |
| def _extend_edges(start: float, limit: float, step: float) -> np.ndarray: |
We never but the step in the middle.
|
|
||
| The focussing step aligns its internal, finer 2theta binning with these edges, | ||
| so setting them also affects results that are integrated over 2theta. | ||
| ``None`` means that the data will not be grouped by 2theta. |
There was a problem hiding this comment.
Nonemeans that the data will not be grouped by 2theta.
That is not true. It means that the pipeline will fail to run if the user requests data that is grouped by two theta. See group_two_theta.
There was a problem hiding this comment.
You are right. Reworded to say that no 2theta grouping is available, i.e. the workflow can only produce results integrated over 2theta and raises if grouped data is requested. The PR description made the same wrong claim and is fixed too.
| KeepEvents[SampleRun](keep_events), | ||
| ) | ||
| grouped = group_two_theta(focussed, two_theta_bins) | ||
| return grouped if grouped.bins is None else grouped.hist() |
There was a problem hiding this comment.
| return grouped if grouped.bins is None else grouped.hist() | |
| return grouped.hist() if grouped.is_binned else grouped |
| [ | ||
| sc.linspace('two_theta', 75.0, 105.0, num=180, unit='deg'), | ||
| sc.linspace('two_theta', 0.8, 2.4, num=17, unit='rad'), | ||
| sc.array(dims=['two_theta'], values=[0.5, 0.6, 1.5, 1.55, 2.9], unit='rad'), |
There was a problem hiding this comment.
What about non-monotonic bins?
- Bins that go up and down: I doubt there is a use case for it. But what does the code do in that case?
- Bins with 0-width: Could happen, does the code cope?
There was a problem hiding this comment.
Good question, both were broken.
Descending edges: focussing died with numpy's ValueError: repeats may not contain negative values. Now the unsorted grid reaches scipp.bin/scipp.hist, which raise BinEdgeError: Bin edges ... must be sorted. I had an explicit check in _focussing_two_theta_bins for a while and then dropped it -- scipp's error is the right type and message, and it fires for both KeepEvents branches. I agree there is no use case, but the failure should at least be legible.
Zero-width bins: an interior one produced a RuntimeWarning from a 0/0 division (which filterwarnings = ["error"] turns into a test failure), and as first or last bin a ValueError: cannot convert float NaN to integer, since the outermost sub-bin width is what continues the grid past the requested range. Such a bin is now subdivided into a single zero-width sub-bin and comes out empty, which is what scipp.hist and groupby do with zero-width bins anyway. Both cases are tested.
There was a problem hiding this comment.
Where are you testing the zero-width interior bin?
There was a problem hiding this comment.
No, I was not. The case had zero-width first and last bins only. Split into zero-width-interior-rad ([0.5, 1.5, 1.5, 2.9]) and zero-width-outermost-rad.
One correction to what I wrote above while I was at it: if I revert the clamp, both cases fail on the 0/0 warning rather than on wrong numbers. A skipped zero-width bin still leaves its edge in the grid as the neighbouring bin's edge, so alignment and counts are unaffected either way, and the NaN no longer reaches the cast because the extension step falls back to the base width when the outermost sub-bin has no width. So what the two cases pin is that no warning or crash is produced, not a numerical difference.
|
I don't quite understand why we have the intermediate theta binning at all? What is the purpose of it? |
Stream processing where we need to be able to accumulate numerator and denominator in histogram-mode. |
Responses to review: - Drop the redundant `max` in `_extend_edges`; `np.arange` already yields no edges when the requested bins extend past the detector range. Order its arguments as start/limit/step. - Subdivide a zero-width requested bin into a single zero-width sub-bin instead of none. It comes out empty, as it would from `scipp.hist`, whereas previously an interior one produced a 0/0 warning and an outermost one a NaN cast error. Descending bins now reach `scipp.bin` and `scipp.hist`, which raise `BinEdgeError`. - Group in the unit of the data and label the result with the requested edges instead of round-tripping the data through the requested unit. The requested edges are then bit-identical to the grid built by focussing, so `group_two_theta` can check exactly that both steps were given the same bins. Without it, calling the two functions directly with different bins silently returns the aliasing this branch removes. - Say accurately what `TwoThetaBins = None` means.
Fair, and the failure was silent, which is the worse part. Making that check exact rather than tolerance-based also removed a rad->deg round trip: grouping runs in the unit of the data now and labels the result with the requested edges, so the output coordinate is the requested one rather than a converted copy of it. |
Aha, it's about the monitor-histpgram normalization. We need the theta-dspacing grid where the monitor correction is applied to be fine enough that each bin corresponds to a small enough range of wavelengths. Is that it? If that's the case I see the problem, but it seems a bit strange to me that the regular non-streaming workflows are this much affected by the requirements of the live-data case. |
I agree that this is not a great situation, but we felt that maintaining two separate workflows would have been even worse. |
The previous case had zero-width first and last bins but no interior one, which is the case that produced the 0/0 warning.
3 participants
Closes #670.
Problem
Focussing binned into 1024 uniform 2theta bins spanning the detector's full range, chosen independently of the 2theta binning requested for the result.
group_two_thetathen assigned each focussing bin as a whole to the requested bin containing its center. Because the two grids are unaligned, the number of focussing bins per requested bin varies periodically, producing large spikes. On a smooth, noise-free 2theta distribution with the binning from the issue (75–105 deg, 180 bins, detector spanning 10–170 deg) the artefact is 93% peak-to-peak;rebinof the same histogram gives 0.05%.This affects the
KeepEvents=Truepath equally:two_thetais a per-pixel coordinate, not an event coordinate, so binning has already discretised it by the time grouping happens.Approach
The focussing grid is now a refinement of
TwoThetaBins: every requested edge occurs exactly in the grid, and each requested bin is subdivided into equal sub-bins until they reach the resolution needed to reconstruct wavelength for monitor normalization (the full 2theta range divided by 1024, i.e. the previous hard-coded value). Grouping into the requested bins is then an exact sum over whole sub-bins — no aliasing, and no fractional splitting of counts asrebinwould introduce.Outside the requested range the grid continues at the base resolution, so
integrate_two_thetastill covers the whole detector and the 1D path keeps its meaning. Subdivision is per requested bin, so a non-uniformTwoThetaBinsdoes not force over-refinement of its narrow bins.rebinwas the cheaper alternative but leaves counts split across unaligned edges, which correlates neighbouring output bins and makes their variances optimistic.Notes
TwoThetaBinsgains aNonedefault, which makes 2theta grouping unavailable: the workflow can then only produce results integrated over 2theta, andgroup_two_thetaraises if grouped data is requested. Workflows that only produce 1D d-spacing output (e.g. the basic DREAM notebook) therefore do not gain a new required parameter.TwoThetaBins, so changing it invalidates accumulated data in a live job. This is unavoidable: an unaligned grid cannot be re-derived exactly. Its size is essentially unchanged for realistic parameters (~1030–1190 bins vs. 1024 before; esslivedata's DREAM default lands around 1100), and only grows if more than ~1024 bins are requested.two_thetacoordinate in radians previously raised, since neithergroupbynorbinconverts units. Focussing now converts the requested bins to the unit of the data, and grouping runs in that unit and labels the result with the requested edges. The output coordinate is therefore identical to what was asked for, rather than a converted round trip of it.group_two_thetachecks that every requested edge occurs in the two-theta grid of the data and raises otherwise. Within the workflow both steps receive the sameTwoThetaBinsand cannot disagree, but a direct caller passing different bins to the two functions would otherwise silently get back the aliasing this PR removes.scipp.histandscipp.groupbydo with such bins anyway. Descending edges are passed on toscipp.bin/scipp.hist, which raiseBinEdgeError.scipp.binandscipp.histkeep their fast path for linear edges. This requires continuing the grid beyond the requested range at the width of the sub-bins rather than at the base width; since a sub-bin is never narrower than half the base width, that at most doubles the number of bins outside the requested range. Losing the fast path cost 12% of the histogramming step here, and would cost a factor of 8-10 if two-theta were an event coordinate rather than a per-pixel one. A non-uniform requested binning cannot produce a linear grid and does not get the fast path.Test plan
New
tests/powder/grouping_test.pycompares the focussed-and-grouped result against a direct histogram of the raw events into the requested bins — the exact reference — for bothKeepEventsbranches and for uniform/non-uniform, degree/radian binning, bins extending beyond the detector range, and zero-width bins.Each behaviour was reverted individually, keeping the new signature so that nothing fails merely on a
TypeError, to check that the tests pin what they claim to:matches_direct_histogramcases,preserves_requested_binsand bothfocussing_with_descending_bins_raisesgroup_two_thetapreserves_requested_binsgroup_two_thetagroup_two_theta_with_misaligned_bins_raiseszero-width-radcases and bothfocussing_with_descending_bins_raisesintegrate_two_theta_covers_full_detector_rangegroup_two_theta_without_requested_bins_raisesfocussing_without_requested_bins_keeps_all_countsNote that the degree binning fails outright, the unit mismatch being caught before grouping, rather than producing wrong numbers. It is a separate defect from the aliasing rather than a symptom of it.
The DREAM and POWGEN reduction suites pass unchanged.