Add waveform arithmetic and statistical reductions - #45
Conversation
Support `+`, `-`, unary `-` between waveforms sharing a time axis, scalar `*`, `/`, `\` on a single waveform, and sample-wise `sum`, `mean`, `var`, `std` across an ArrayOfRDWaveforms. Reductions require all waveforms to share a time axis and throw an ArgumentError otherwise; the shared axis carries through to the result. Makes superpulse construction a direct `mean(wfs)` rather than requiring a dedicated function. Assisted-by: Claude Opus 5 <noreply@anthropic.com>
Sample-wise reductions accumulate into a single preallocated buffer instead of combining whole signal vectors pairwise, which makes their cost independent of how signals are stored: unchanged for contiguous storage, and roughly twelve times faster with a fraction of the memory for signals held as separately allocated vectors. Integer samples narrower than Int now accumulate in Int, since summing thousands of Int32 detector samples overflows silently. Broadcasting an operator over an ArrayOfRDWaveforms is evaluated eagerly over the whole sample storage, so contiguously stored signals stay contiguous rather than being rebuilt as a vector of separately allocated vectors. This trades away broadcast fusion, and matches how RadiationDetectorDSP treats the equivalent named operations. Assisted-by: Claude Opus 5 <noreply@anthropic.com>
Shifting every sample by a constant is one of the most common waveform operations, but was reachable only by constructing a second waveform. Support `wf + a` and `wf - a` on a single waveform, and broadcast forms on an ArrayOfRDWaveforms taking either one shift for all waveforms or a vector holding one shift per waveform, which covers subtracting per-waveform baselines. The scalar is typed as a RealQuantity, so adding a bare number to a waveform carrying units is rejected rather than silently accepted. `wfs + a` on an ArrayOfRDWaveforms is intentionally left undefined: it is an AbstractArray, and adding a scalar to an array without broadcasting is not valid Julia. Assisted-by: Claude Opus 5 <noreply@anthropic.com>
Shifting a waveform whose samples carry a unit by a plain number is unambiguous: the number is the shift expressed in that same unit. It now does that instead of being rejected as a dimension error, which matters because sample units are often attached only partway through an analysis while the shifts computed from the data stay plain numbers. Scalars that carry their own unit keep converting as before, and waveforms with unitless samples are unaffected. Assisted-by: Claude Opus 5 <noreply@anthropic.com>
Shifting combines samples and a shift amount, and either may carry a unit while the other does not. In that case the side without one is now taken to be expressed in the unit of the other, in both directions, so plain samples accept a unitful shift just as unitful samples accept a plain one. Sample units are often attached only partway through an analysis, and neither mixture is ambiguous. Shifts that already agree in unit convert as before, and shifting unitless samples by a plain number leaves them unitless. Assisted-by: Claude Opus 5 <noreply@anthropic.com>
The sample-wise mean, var and std add methods to the generic functions Statistics defines, so name Statistics as the dependency they come from and declare it in Project.toml. StatsBase re-exports those same generics, so the methods themselves are unchanged. Assisted-by: Claude Opus 5 <noreply@anthropic.com>
The suite had grown to cover the same behaviour several times over: scaling and division were retested verbatim with units even though they do not touch units, every broadcast operator ran against both signal storage layouts although the layout only selects one shared dispatch path, and the unit-matching rules were spread over three testsets. Group the broadcast forms into a table checked in one loop, exercise the generic storage path with one operator of each kind, and state each unit-matching rule once. Same behaviours covered, a third fewer checks. Assisted-by: Claude Opus 5 <noreply@anthropic.com>
Sample-wise sum/mean/var/std now specialize on contiguous storage, reducing over flatview(signals) in one block instead of one broadcast per waveform. This is a modest win on a CPU and the deciding factor on a device that dispatches each pass separately, where the per-waveform loop runs tens of times slower.
JLArrays.jl provides an AbstractArray backend with no special casing anywhere in Base or this package, the same property a real GPU array type (CuArray, ROCArray, ...) has. Add it as a test dependency and exercise the ArrayOfSimilarVectors/flatview-specialized reductions and every broadcast form against it, confirming the fast path dispatches correctly and preserves non-Array storage without requiring GPU hardware.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #45 +/- ##
===========================================
+ Coverage 12.38% 54.41% +42.02%
===========================================
Files 4 4
Lines 113 215 +102
===========================================
+ Hits 14 117 +103
+ Misses 99 98 -1 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Statistics only gained registry-tracked versions at 1.11.0, so a "1" compat bound has no satisfiable floor and fails downgrade CI, which resolves every direct dependency to its lowest compatible version. Statistics is a stdlib bundled with every Julia version this package supports (down to the 1.10 LTS floor); leaving it uncompatted lets it resolve to whatever ships with the running Julia, same as the package's other stdlib dependencies.
|
Failing doc build might be fixed by #46. |
|
ChatGPT 5.6 Sol max review (@theHenks you can forward this to your coding AI): Review verdict: request changesThe PR is generally well put together: the API additions are coherent, time-axis compatibility is checked, the reduction implementation pays attention to integer widening, contiguous storage, nonstandard axes, units, and non-
Aside from those points, I like the implementation direction. The two-pass variance is preferable to So my merge call would be: fix/test heterogeneous ragged reductions, then approve. I would treat the eager-broadcast and unit semantics as follow-up/documentation items rather than merge blockers. I reviewed the current PR branch/source and tests statically; I did not execute the Julia test suite locally because Julia isn't available in this environment. |
docs/Project.toml pinned Documenter to "~0.24" (resolving to 0.24.11) while test/Project.toml already used "~0.27" for the doctest step. Documenter 0.24's docstring-coverage check cannot correctly attribute docstrings on methods that extend an external function (Statistics.mean/ var/std here) with a parametric where-clause signature, so it flagged them as missing and failed the strict docs build. Verified locally: with Documenter bumped to 0.27.25, the full strict build (CheckDocument through HTMLWriter) completes with no missing-docstring warning. Also declares RadiationDetectorSignals itself as a docs dependency (Pkg.develop path-installed), matching what the CI docs job already does before running make.jl.
_sample_sum derived its accumulator's element type from only the first signal, so a ragged collection mixing e.g. an Int32 waveform with a Float64 one would size the buffer to Int32/Int and lose or fail to represent later, wider-typed contributions. Derive it from the promoted type across every signal instead, matching Julia's normal mixed-argument promotion. sum/mean/var/std on ArrayOfRDWaveforms all inherit the fix, since they build on _sample_sum. The flatview- specialized path for contiguous storage is unaffected: a single contiguous block can only ever hold one concrete element type. Reported by @oschulz.
The +/- scalar-shift docstrings on RDWaveform only documented the "plain number takes the samples' unit" direction. State both directions explicitly: whichever of the samples and the shift amount carries no unit takes on the other's. Requested by @oschulz.
This reverts commit 70b1594.
Statistics only gained registry-tracked versions at 1.11.0, so downgrading its compat entry fails no matter which Julia version runs the job: the action narrows "1" to "~1.0.0", and neither the registry range (1.11.0-1.11.1, checked on Julia "1"/latest) nor the bundled stdlib version (1.10.0 on Julia 1.10) falls inside that window. Verified both independently by reproducing the action's exact rewrite and resolving/building against it on each Julia version. Running the job on 1.10 instead of "1" also tests against the package's declared LTS floor (`julia = "1.10"`) rather than whatever Julia currently calls latest. Skipping Statistics from the downgrade leaves its compat at "1", which resolves fine on 1.10 (bundled Statistics reports 1.10.0, well within [1.0.0, 2.0.0)). Verified end-to-end: reproduced the action's actual compat rewrite for every dependency, resolved and built against it on Julia 1.10, and ran the full test suite — 159/159 passing.
A unitful shift amount on a plain waveform (wf + 10u"ns") silently gave the samples that unit rather than shifting them, which reads as a time shift but actually reinterprets an ADC/energy signal as being measured in nanoseconds. Drop that direction: a plain shift is still interpreted in unitful samples' own unit, but a unitful shift on plain samples is now a Unitful.DimensionError, the same as adding any two incompatible Unitful quantities, rather than a silent unit promotion. Flagged by @fhagemann.
The reduction docstrings only mentioned the ArgumentError thrown for mismatched time axes, not that a DimensionMismatch is also possible when the signal arrays themselves don't share axes (a different check, inside the sample-wise accumulation). Flagged by @fhagemann.
The offset-axes generic-indexing coverage rolled its own ~20-line ShiftedVector/BroadcastStyle pair to get a vector with non-1-based indices. OffsetArrays.jl is the standard package for exactly this and is already well-tested; use its OffsetVector as a test-only dependency instead of duplicating what it does. Flagged by @fhagemann.
Nothing in the package uses StatsBase: the only reference was a Histogram method left commented out at the end of detector_waveforms.jl, predating this branch. The commented sketch stays as a placeholder for whoever picks that up; only the live import and the compat/deps entries are removed.
|
@oschulz @fhagemann This is ready from my side. |
The arithmetic and shift operators carry no docstrings of their own; the capabilities they provide are documented once, on the type.
Shifting a waveform by a scalar uses plain broadcasting: the amount must be dimensionally compatible with the samples, and mixing a unitful amount with plain samples (or the reverse) throws a Unitful.DimensionError.
`Base.promote_op(Base.add_sum, T, T)` is the element type `sum` itself would produce, so narrow-integer widening no longer depends on an enumerated list of sample types.
Equal-length signals held in a VectorOfVectors reshape into a matrix without copying, so they reduce in one pass over their own storage like an ArrayOfSimilarVectors already did. Only signals with no contiguous block behind them fall back to accumulating waveform by waveform. Variance on that fallback defers to Statistics.varm, given the sample-wise mean so that narrow integer samples still accumulate without overflow.
Broadcasting and per-waveform shifts take the same matrix view of the samples that the reductions do, so equal-length signals held in a VectorOfVectors stay in one block instead of being rebuilt as separately allocated vectors.
Each operation in a composed expression materializes in turn, so expressions combining waveform, scalar and per-waveform operands are checked against the same computation applied waveform by waveform, across every signal storage layout and on the JLArrays backend.
Int32 is narrower than Int only where Int is 64 bits, so the widening is checked with Int16 on every platform and with Int32 where it applies.
`sum`, `mean`, `var` and `std` are heavily specialized across packages, so their behavior for waveform arrays is described on the type rather than attached to the functions themselves. The same docstring now also covers broadcast arithmetic over an array of waveforms. Assisted-by: Claude Opus 5 <noreply@anthropic.com>
| Base.float(wf::RDWaveform) = RDWaveform(float(wf.time), float(wf.signal)) | ||
|
|
||
| function Base.:(+)(a::RDWaveform, b::RDWaveform) | ||
| a.time == b.time || throw(ArgumentError("Can't add RDWaveform with different time axes")) |
There was a problem hiding this comment.
Maybe an approximate comparison instead of == ?
| end | ||
|
|
||
| function Base.:(-)(a::RDWaveform, b::RDWaveform) | ||
| a.time == b.time || throw(ArgumentError("Can't subtract RDWaveform with different time axes")) |
There was a problem hiding this comment.
Maybe an approximate comparison instead of ==?
|
|
||
| Base.:(-)(a::RDWaveform) = RDWaveform(a.time, -a.signal) | ||
|
|
||
| Base.:(+)(wf::RDWaveform, a::RealQuantity) = RDWaveform(wf.time, wf.signal .+ a) |
There was a problem hiding this comment.
Technically this should probably be a broadcasted plus. But making semantically clean and work robustly via broadcast specialization would be non-trivial (see comment on broadcasted specialization s below). S I think going for Base.:(+) could be tolerable - it's mathematically a bit ugly, though:
If we define algebra on waveforms, we should define what waveforms are, mathematically, with respect to algebra operations. We could see waveforms with equal time axis as elements of a vector space: then addition of waveforms and multiplication of waveforms with scalars would be well defined, and so we be operations like sum, mean, etc. But addition between waveforms and scalar values would not be well defined, at least not without additinal mathematical gymnastics/definitions.
There is a clean solution though: We introduce a special signal type RDBaseline(level) that represents an axis-independent constant signal with amplitude level. Then wf + RDBaseline(level) would be both mathematically well-defined and also very readable without digging into the docs of RadiationDetectorSignals.
We would define
"""
RDBaseline(level)
An axis-independent constant signal with amplitude `level`.
When combined with an `RDWaveform`, it represents the constant waveform with
value `level` on that waveform's time axis.
"""
struct RDBaseline{T<:RealQuantity}
level::T
end
Base.:(+)(bl::RDBaseline, a::RealQuantity) = RDBaseline(bl.level + a)
Base.:(+)(a::RealQuantity, bl::RDBaseline) = bl + a
Base.:(-)(bl::RDBaseline, a::RealQuantity) = RDBaseline(bl.level - a)
Base.:(-)(a::RealQuantity, bl::RDBaseline) = RDBaseline(a - bl.level)
Base.:(*)(a::Real, bl::RDBaseline) = RDBaseline(a * bl.level)
Base.:(*)(bl::RDBaseline, a::Real) = a * bl
Base.:(/)(bl::RDBaseline, a::Real) = bl * inv(a)
Base.:(\)(a::Real, bl::RDBaseline) = bl / a
Base.:(+)(wf::RDWaveform, bl::RDBaseline) = RDWaveform(wf.time, wf.signal .+ bl.level)
Base.:(+)(bl::RDBaseline, wf::RDWaveform) = wf + bl
Base.:(-)(wf::RDWaveform, bl::RDBaseline) = RDWaveform(wf.time, wf.signal .- bl.level)
Base.:(-)(bl::RDBaseline, wf::RDWaveform) = RDWaveform(wf.time, bl.level .- wf.signal)(untested, may contain bugs).
A clear distinction between a generic scalar quantity and a constant waveform could be very useful in general.
| ArrayOfRDWaveforms((a.time, _broadcast_signals(f, a.signal, b.signal))) | ||
| end | ||
|
|
||
| Base.Broadcast.broadcasted(::typeof(*), a::Real, wfs::ArrayOfRDWaveforms) = |
There was a problem hiding this comment.
Is this robust? What if wfs is a Broadcasted, due to broadcast-chaining - we'd have different semantics or an error then, right?
Maybe not do this for now? It would be weird to define a .* wf if a * wf[i] is not well defined. A robust solution would require an iteration scheme over waveforms as custom pairs or time and signal values, with scalar multiplication defined for the pair (action only on the signal value) and a specialized broadcast output scheme.
| return out | ||
| end | ||
|
|
||
| function _sample_sum(signals::AbstractVector{<:AbstractVector}) |
There was a problem hiding this comment.
We should add something like
_sample_sum(signals::ArrayOfSimilarVectors) = sum(signals)
since ArraysOfArrays provides optimized and GPU-compatible specializations for sum,mean, var, etc.
| # storage; `nothing` when the signals are not one contiguous block of equal-length | ||
| # vectors. Operating on the whole block in a single pass keeps the work on whatever | ||
| # device holds the samples, instead of dispatching one pass per waveform. | ||
| _sample_matrix(signals::ArrayOfSimilarVectors) = flatview(signals) |
There was a problem hiding this comment.
If we require ArraysOfArrays v1, we can just use stacked(signals) instead of _sample_matrix.
| Base.Broadcast.broadcasted(::typeof(-), wfs::ArrayOfRDWaveforms) = | ||
| _scaled_waveforms(wfs, x -> .-x) | ||
|
|
||
| Base.Broadcast.broadcasted(::typeof(+), a::ArrayOfRDWaveforms, b::ArrayOfRDWaveforms) = |
There was a problem hiding this comment.
Why a different implementation than Base.:(+)(a::RDWaveform, b::RDWaveform)?
Also, this may not be a good idea, see comment on Base.Broadcast.broadcasted(::typeof(*), ...).
| function _nested_sample_sum(signals::AbstractVector{<:AbstractVector}) | ||
| out = similar(first(signals), _sample_sum_eltype(eltype(eltype(signals)))) | ||
| fill!(out, zero(eltype(out))) | ||
| for signal in signals |
There was a problem hiding this comment.
Summation via an accumulator can be numerically unstable, especially when summing over a large number of entries. sum(signals) uses advanced summation schemes and is numerically stable, but generates a lot of memory allocations.
A good approach might be to use sum(nestedview(stack(signals))) (or better sum(sliced(stacked(signals))) which would just do (at most) one or two memory block allocations/copies and use numerically stable summation.
We could define a helper function
_aosa(vs::VectorOfSimilarVectors} = vs
_aosa(vs::AbstractVector{<:AbstractVector}} = sliced(stacked(signals)) # nestedview(stack(signals)) for ArraysOfArrays v0.6and then use that in waveform sums, means, etc. and let ArraysOfArrays use its specializations.
|
|
||
| `sum`, `mean`, `var` and `std` reduce sample-wise over all waveforms and | ||
| return a single [`RDWaveform`](@ref) carrying the shared time axis. `var` and | ||
| `std` use the Bessel-corrected denominator `length(wfs) - 1`. Integer samples |
There was a problem hiding this comment.
| `std` use the Bessel-corrected denominator `length(wfs) - 1`. Integer samples | |
| `std` resprent the pointwise variance and standard error and use the Bessel-corrected denominator `length(wfs) - 1`. Integer samples |
3 participants
Disclaimer: I applied my newly gained claude code knowledge from the juliacon to this PR!
This resurrects and extends the stashed
waveform-mathbranch, which came up in the discussion on JuliaPhysics/RadiationDetectorDSP.jl#71 (sum_waveforms) as the more natural home for this kind of operation.What's added
RDWaveform:+,-, unary-between waveforms (matching time axes required); scalar*,/,\; scalar+/-shifts in both argument orders, with either side free to carry a unitArrayOfRDWaveforms:sum,mean,var,stdsample-wise across waveforms, including aflatview-based fast path for contiguous storageBase.Broadcast.broadcastedspecializations so broadcasting the above over anArrayOfRDWaveformspreserves the underlying signal storage layout (contiguous storage stays contiguous), following the same patternRadiationDetectorDSP.jluses forshift_waveform/multiply_waveformTesting
Pkg.test()passes 153/153, including the Documenter doctestsJLArrays.jl(a non-ArrayAbstractArraybackend used for GPU-generic testing without requiring GPU hardware), confirming the storage-preserving fast paths dispatch and produce correct results on a backend other thanArraymean/sum)Compatibility
Statisticsis added as a new direct dependency (stdlib, available on the 1.10 LTS floor) for themean/var/stdextensions. No breaking changes.CC @oschulz — saw your note on #71 about possibly tackling this yourself; just wanted to have a look and test my claude code setup. let me know if you intended something else here