You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
{{ message }}
Repository navigation
perf(metrics): fused GPU kernels and batched MicroMS3IM, GLCM and spectral sums - #58
GPU speedups for the cubic metrics the dynacell eval spends its non-segmentation time in. There are two
kinds of change:
fused CUDA reductions that replace chains of full-size array operations or contended bincount atomics;
batched entry points that replace per-slice and per-cell Python loops.
Most values are unchanged or move at 1e-16. The one larger change is Spectral_PCC, at 4.3e-7 relative
(details below). The PR adds two new APIs and changes no existing signature.
Motivation: per FOV, the A549 ER eval on an A40 spends 7.2 s in MicroMS3IM scoring, 7.1 s in its
calibration, about 6 s in per-cell GLCM and 4.7 s in pixel metrics (dynacell eval_timing.csv and py-spy).
These follow the GPU segmentation in #57.
FRC/FSC and radial-power bin sums in one shared-memory kernel instead of 3–4 weighted bincounts (opt-in shared memory; out-of-range ids, host ids and sm < 60 take the bincount path)
dynacell compute_pixel_metrics: 403.6 → 299.5 ms per timepoint
The fit and the scoring together, against origin/main on the same inputs: fit 45.9 → 1.48 s (31×), α
identical (3.2737876176834106).
Numerical equivalence
MicroMS3IM:
α is identical.
Per-timepoint scores are identical to main and to the production run's.
Per-slice scores differ by ≤ 2.2e-16.
Per-label GLCM:
every property is within 3.4e-16 relative of the per-cell glcm_features loop;
per-region ranges work in float64, so integer, bool and NaN-containing images match glcm_features too;
sparse and negative label ids are supported.
Pixel metrics:
PCC, SSIM, NRMSE, PSNR and their SI variants, and both FSC resolutions, are identical.
Spectral_PCC moves by 4.3e-7 relative (7.6e-9 absolute), because reduce_power on GPU forms |F|²
as re² + im² rather than squaring |F|.
The fused kernels compute each per-element term with the array path's float operations in the same order,
compiled with --fmad=false. Only the float64 summation order differs. The bincount path they replace
also added atomically, in an arbitrary order.
Tests
Fused RI objective:
matches the CuPy maps for float32 and float64, in 3D, 2D and 1D layouts;
a spy asserts that the fused path ran;
the GPU fit matches the CPU fit.
Fixed-point exit:
α equals the old loop run to its cap;
the test asserts that the old loop never converged on that input.
score_stack:
equals per-slice score on CPU and GPU, across three batch sizes;
includes constant-GT slices and a batch with no valid slice;
covers input validation.
glcm_features_by_label:
equals glcm_features with each region's mask, in 2D/3D, on CPU/GPU, with shared and per-region
ranges;
covers adjacent labels, a hole, a constant region and a 1-voxel region (NaN);
a mutation that counts cross-label pairs fails it.
Fused spectral sums:
match the bincount path for complex64 and complex128, plain, sectioned and power-only;
unsupported inputs are declined.
Test plan
ruff check ., ruff format --check ., mypy --ignore-missing-imports cubic/
full pytest on an A40: 862 passed, 2 skipped (exit 0)
/simplify, Copilot, /review and /code-review findings addressed: one commit per finding, threads resolved
Accelerate GPU metric evaluation with fused reductions and batched scoring while adding batched MicroMS3IM and per-label GLCM APIs.
New Features:
Add batched MicroMS3IM stack scoring with per-slice handling of degenerate inputs.
Add batched per-label GLCM feature extraction for sparse, negative, and mixed-size labels.
Bug Fixes:
Stop RI-factor bisection when its bounds reach a floating-point fixed point.
Preserve correct spectral reductions by falling back for unsupported or out-of-range GPU inputs.
Enhancements:
Fuse GPU reductions for RI-factor objectives and spectral power/cross-spectrum sums to accelerate metric computation.
Refactor MS-SSIM internals to support per-image data ranges for batched scoring.
Extend GLCM quantization and validation to support per-region ranges and shared NumPy/CuPy execution.
Documentation:
Document the new batched MicroMS3IM and per-label GLCM APIs and their input behavior.
Tests:
Add CPU and GPU coverage for fused RI-factor reductions, batched MS-SSIM scoring, per-label GLCM extraction, and fused spectral sums, including validation and fallback cases.
Each mean S / dS evaluation built ~15 element-sized temporaries over the
pooled elements. On GPU, float32/float64 elements now go through a CuPy
ReductionKernel that repeats the NumPy operation order with FMA
contraction off, so the per-pixel values equal the array path's and only
the float64 summation order changes. Rows of the leading axis reduce in
separate blocks, since a reduction to one scalar runs on a single block.
MEASURED on an A40, 576 captured 634x954 float32 slices: one dS
evaluation 221 ms -> 15.6 ms, alpha unchanged.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
The loop stops only when |f| < 1e-10 and |hi - lo| < 1e-8, but with
float32 elements |f| never drops below 1e-10. The loop then ran to
_MAX_BISECT_ITERS long after lo and hi became adjacent floats, where
every iteration repeats the same state. Break there instead: the
returned alpha is identical, and the captured 576-slice fit pool
needs 56 objective evaluations instead of 204.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
…M.score_stack
Scoring an FOV called MicroMS3IM.score once per z-slice: ~15 small
filter and reduction launches per scale, and a host sync for every
slice's data range. score_stack normalizes and scores up to 2**24 pixels
of slices per pass, each slice with its own ground-truth data range
(ms_ssim's per-image path now accepts an (N,) data_range). Slices whose
range is not finite and positive, where score raises, get a caller-chosen
fill value.
MEASURED on an A40, one captured FOV (480 slices of 640x960), against the
per-slice loop on origin/main: 4.92 s -> 0.72 s. Every slice agrees to
2.2e-16, and the per-timepoint means are identical.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
glcm_features takes one region per call, so per-cell texture meant one
call per cell: a bounding-box crop, 13 direction passes and a host sync
for each. glcm_features_by_label accumulates the co-occurrences of all
regions in one bincount per direction over the whole image, keyed by
(label, direction, i, j), and computes the Haralick properties vectorized
over labels. Each region equals glcm_features with that region's mask;
a region without any voxel pair is NaN instead of raising.
MEASURED on an A40, one captured A549 ER timepoint (48x640x960, 13
cells, levels=32): dynacell's per-cell loop 0.604 s -> 0.131 s, every
property within 3.4e-16 relative.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
Each weighted bincount in the FRC/FSC and radial-power reducers made
every voxel add atomically in float64 into a handful of global bins.
Those adds contended: one bincount took 32 ms on a 48x640x960 spectrum,
and a sectioned FSC needs three. Complex GPU spectra now go through one
kernel. It computes |X|^2, |Y|^2 and Re{X conj(Y)} per voxel with the
array path's float ops, accumulates them per block in shared memory,
and flushes once. Binnings that do not fit in 48 KB of shared memory,
host arrays and other dtypes keep the bincount path.
MEASURED on an A40, one captured A549 timepoint through dynacell
compute_pixel_metrics: 403.6 -> 298.1 ms. Both FSC resolutions are
identical; Spectral_PCC moves by 4.3e-7 relative because reduce_power
now forms |F|^2 as re^2 + im^2 instead of squaring |F|.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
The PR adds two public batched APIs and GPU-specific fused reductions to remove full-array intermediates, contended histogram atomics, and per-slice/per-label Python loops. It preserves existing signatures through guarded fast paths and fallbacks, with targeted tests covering numerical equivalence, validation, degenerate inputs, GPU dispatch, and unsupported-kernel cases.
Trigger a new review: Comment @sourcery-ai review on the pull request.
Continue discussions: Reply directly to Sourcery's review comments.
Generate a GitHub issue from a review comment: Ask Sourcery to create an
issue from a review comment by replying to it. You can also reply to a
review comment with @sourcery-ai issue to create an issue from it.
Generate a pull request title: Write @sourcery-ai anywhere in the pull
request title to generate a title at any time. You can also comment @sourcery-ai title on the pull request to (re-)generate the title at any time.
Generate a pull request summary: Write @sourcery-ai summary anywhere in
the pull request body to generate a PR summary at any time exactly where you
want it. You can also comment @sourcery-ai summary on the pull request to
(re-)generate the summary at any time.
Generate reviewer's guide: Comment @sourcery-ai guide on the pull
request to (re-)generate the reviewer's guide at any time.
Resolve all Sourcery comments: Comment @sourcery-ai resolve on the
pull request to resolve all Sourcery comments. Useful if you've already
addressed all the comments and don't want to see them anymore.
Dismiss all Sourcery reviews: Comment @sourcery-ai dismiss on the pull
request to dismiss all existing Sourcery reviews. Especially useful if you
want to start fresh with a new review - don't forget to comment @sourcery-ai review to trigger a new review!
The reason will be displayed to describe this comment to others. Learn more.
Hey - I've found 2 issues
Prompt for AI Agents
Please address the comments from this code review:
## Individual Comments### Comment 1
<locationpath="cubic/metrics/microssim/micro_ms3im.py"line_range="154-156" />
<code_context>
+ K2=K2,+ normalize=normalize,+ )+ scores[start : start + len(valid)][valid] = asnumpy(+ cast(np.ndarray, batch_scores)+ )+ return scores+
</code_context>
<issue_to_address>
**issue (bug_risk):**`scores[start : start + len(valid)][valid] = ...` assigns through a boolean-indexed copy, so the original `scores` array is never updated. Every valid slice therefore retains the initial `degenerate` value instead of receiving its MS-SSIM score.
**Triggers:** Whenever `score_stack` contains at least one valid slice.
**Suggested fix:** Assign through a single indexed expression, such as `scores[start : start + len(valid)][valid]` via an explicit temporary index (`scores[start + np.flatnonzero(valid)] = ...`) or assign the filtered result back to the slice.
```suggestion scores[start + np.flatnonzero(valid)] = asnumpy( cast(np.ndarray, batch_scores) )```
</issue_to_address>
### Comment 2
<locationpath="cubic/feature/texture.py"line_range="443-445" />
<code_context>
+# Compact labels 1..n_labels in ascending label order; 0 stays background.+ xp = get_array_module(labels)
+ index = xp.arange(1, n_labels + 1)
+ lookup = xp.zeros(int(label_ids[-1]) + 1, dtype=np.int64)
+ lookup[present] = index
+ compact = lookup[labels]+ del lookup
+
</code_context>
<issue_to_address>
**issue (performance):** The compact-label lookup allocates an array of length `max(label) + 1`, so a valid sparse label image such as one containing label `2**31` attempts a multi-gigabyte allocation before any GLCM work and can raise an allocation error or exhaust GPU memory.
**Triggers:** When labels are sparse or contain a large integer identifier.
**Suggested fix:** Compact labels using a sort/search mapping or a dictionary-like indexed operation that scales with the number of distinct labels rather than the maximum label value.
</issue_to_address>
_chunked_mean picked the fused kernel by looking its map function up in
a table, so the generic chunked reducer carried a special case and a
wrapped map silently lost the fast path. _compute_S_mean and
_compute_dS_mean now dispatch to the kernel by name, and _chunked_mean
is the plain array path again. The fused reduction also runs over the
trailing axes instead of reshaping, which copied cropped (strided)
element views on every evaluation.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
…red memory
The kernel was held to the 48 KB default, about 1536 bins, so finer
sectioned FSC binnings (e.g. 480 radial x 6 angular bins without spacing)
fell back to the contended bincount path. Enable each device's opt-in
maximum on the kernel and size the fallback check from it. Also derive
the bin count from n_radial x n_angle instead of passing both, keep the
three sums in one zeroed buffer, and return the counts as an int64 view
rather than a converted copy.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
… entry points
glcm_features_by_label repeated _quantize's floor/clip rule inline, as
well as glcm_features' argument checks and offset list. _quantize now
accepts per-voxel bounds (scalar bounds take the unchanged path), and
both functions call _validate_glcm_args and _glcm_offsets, so
glcm_features' output is unchanged. Per-region extrema use
np.minimum.at / np.maximum.at, which dispatch to CuPy: 28 ms for both
against 445 ms for one ndimage.minimum (48x640x960, 300 labels, A40).
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
… does
score_stack spelled out ms_ssim's defaults (kernel_size=11, sigma=1.5,
...) and copied _prepare's fitted/shape checks, so a change to either
would make it disagree with score. It now takes **ms_ssim_kwargs merged
over defaults read from ms_ssim's signature, rejects unknown names and
data_range, and shares MicroSSIM._check_pair with _prepare.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
…ookup
glcm_features_by_label mapped labels through an array of length
max(label) + 1, so a sparse id such as 2**40 tried a terabyte-scale
allocation and a negative id indexed from the end and collided with a
positive label. Search the sorted present labels instead, which keeps
memory proportional to the image. A test with sparse and negative ids
checks that every region keeps its features.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
The kernel reads both spectra up to FX.size, so a same-dtype but shorter
FY took the fused path and was read past its end. Return None for it, so
the array path raises as before.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
…y 6.0
The fused bin-sum kernel adds doubles atomically in shared memory, which
needs sm_60; on older devices it failed to compile instead of falling
back. Check the device's compute capability before using it.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
glcm_features_by_label combined the labels' device arrays with the
image's in every direction, so a host label image with a GPU image
failed later with an opaque CuPy TypeError. Check both are on the same
device up front.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
The per-region (min, max) kept the image dtype and started from the
whole image's extremes. As a result:
- int16 ranges wider than 32767 overflowed, sending every voxel to level 0;
- bool images raised on subtraction;
- CuPy's ufunc.at rejected uint8/uint16/int16 on the GPU;
- a single NaN anywhere, even in the background, poisoned every region.
Scatter the extremes in float64 from +inf/-inf, as glcm_features works
in Python floats, so a NaN only affects its own region.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
…e properties
Each direction's bincount spanned the bins of all labels and all
directions, so with 13 directions it allocated and added a histogram 13x
larger than the direction needed. The float64 property stage also held
~10 temporaries of the full (labels, directions, L, L) size. Count each
direction into its own slice, and compute the properties over bounded
label chunks. The docstring no longer claims the cost is independent
of the number of regions.
MEASURED on an A40, one captured A549 timepoint (13 cells): 0.134 ->
0.111 s, values unchanged.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
The buffer was declared as unsigned char and reinterpreted as double and
uint64 arrays, so the atomics' 8-byte alignment rested on where the
toolchain happens to place dynamic shared memory. Declare it
__align__(8).
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
…binning
The fused kernel skipped only negative ids, so an id at or above nbins
(or a radial id at or above n_radial in the sectioned form) was added
into the neighbouring sum's shared-memory slots. The array path's
bincount extends its output instead. MEASURED: reduce_power with ids up
to 11 and nbins=10 counted 813 of 1000 voxels into a (10,) result, while
the array path returns (12,) with all 1000. The kernel now counts such
ids, and the caller takes the array path when any occur.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
With a GPU spectrum and host bin ids, the fused path copied the ids to
the device on every call, while the array path fails on the mixed
devices. Whether such a call worked therefore depended on the dtype and
the binning size. Take the array path unless the ids are on the GPU.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
Every _fused_frc_sums call rebuilt CuPy's full device-attribute dict to
read the multiprocessor count, and queried the compute capability again.
Both are fixed per device, so the cached kernel builder now returns them
(and None below sm_60) alongside the kernel and its shared-memory limit.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
score_stack called the per-image MS-SSIM directly and so skipped the
check_same_device guard ms_ssim applies. A GPU gt with a host pred
therefore failed deep in CuPy with a TypeError, where score raises the
library's ValueError.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
_fused_eligible looked only at ux's device and dtype, but the reduction
kernel types all five element arrays alike. Elements of mixed float
dtypes, which the array path handles, therefore raised a CuPy type error.
Require every array to be on the GPU with ux's dtype.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
Only the property stage was chunked: the device histogram and its host
copy still held every label's (direction, i, j) counts at once. With the
3-D defaults that is ~9.9 GiB per copy for 100,000 labels, even when
they are singletons. Count, transfer and reduce each chunk of labels in
turn, holding at most 2**23 histogram bins. A typical FOV is one chunk;
the captured 13-cell timepoint takes 0.123 s (was 0.111 s), values
unchanged.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
a * n_radial + r was computed in int, so a large angle id could wrap to
a negative bin, slip past the bounds check and index shared memory out
of range. Compute it as long long, so such ids take the bincount
fallback.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
This fails on an otherwise supported pre-Pascal GPU: _frc_sums_kernel intentionally declines compute capability < 6.0, so all three GPU calls take the fallback and sum(calls) is zero. Gate the fused-path assertion on the kernel being supported by the current device.
This issue also appears on line 261 of the same file.
Only FX's device was checked. A host FY was copied silently, and an
array on another GPU would have been read by a kernel launched on FX's
device. Take the array path unless FY and the bin ids are CuPy arrays on
FX's device.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
…bility 6.0
On those devices the fused bin sums deliberately take the bincount path.
The tests that assert the fused path ran, or read its shared-memory
limit, would therefore fail there. Skip them through one helper.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
Copilot's "previously missed" note (tests/metrics/spectral/test_radial.py:246, gate the fused-path assertion on GPU kernel support) is fixed in 04fcb57. A _require_fused_kernel() helper skips the tests that assert the fused path ran or read its shared-memory limit on GPUs below compute capability 6.0, where the code takes the bincount path on purpose. Full suite on an A40: 863 passed, 2 skipped.
Reviewed the entire PR diff against its base on main (50b1cbe1efd7df9c3e1c049fe9bd7a1e496c0a08), through head commit 04fcb57ce0d5b854d2b0d3b4d32b2a0cdb7864ce.
MEASURED: With the default per-region range, inserting a NaN into a varied region makes glcm_features_by_label return the exact constant-region feature row on CPU: contrast 0, ASM 1, entropy 0. The same GPU input raises ValueError: The first argument of bincount must be non-negative; the existing per-region glcm_features also raises. This violates the documented equivalence and produces valid-looking texture values for an invalid region.
MEASURED controls: Before inserting the NaN, the affected region's contrast is 0.75 on both devices. An intentionally constant region has contrast 0; the corrupted CPU row matches it for every property. The unaffected region retains its finite-control values. The reproduction and control assertions exited 0 against the PR checkout.
Please reject non-finite labeled regions consistently, or return NaN feature rows for those labels, while continuing to allow NaNs in background voxels. Add coverage for a NaN inside a label on CPU and GPU; the current test covers background NaNs.
MEASURED validation:/hpc/mydata/alex.kalinin/cubic/.venv/bin/python -m pytest -q from the PR checkout on an NVIDIA A40: 863 passed, 2 skipped, 35 warnings, exit 0. gh pr checks 58 --repo alxndrkalinin/cubic also exited 0 with CI passing. No other actionable findings were confirmed.
With per-region ranges, a NaN inside a region made its range NaN. On
CPU the region quantized to one level and returned the constant-region
row (contrast 0, ASM 1, entropy 0), a valid-looking value for an invalid
region; on GPU bincount raised on the cast NaN. Treat a region
containing NaN or inf as background, for per-region and shared ranges
alike, so it has no pairs and its row is NaN. Non-finite background
voxels remain ignored.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
Fix:glcm_features_by_label now returns a NaN feature row for any region that contains a NaN or inf voxel. This holds for both per-region and shared value_range, on CPU and GPU. The region is counted as background, so it has no voxel pairs, which is already the NaN case. Non-finite background voxels are still ignored, and the docstring states this contract.
Test:test_glcm_by_label_nonfinite_region_is_nan puts a NaN or inf inside label 3, for both range modes, on CPU and GPU.
Label 3's row must be all NaN.
The other regions must match glcm_features with their own masks (rel 1e-12).
Against the pre-fix texture.py the test fails. On CPU that code returned the constant-region row; on GPU it raised in bincount.
Verification on an A40:
Full pytest: 871 passed, 2 skipped.
The captured A549 timepoint: 0.124 s (was 0.123 s), values unchanged.
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
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.
Summary
GPU speedups for the cubic metrics the dynacell eval spends its non-segmentation time in. There are two
kinds of change:
bincountatomics;Most values are unchanged or move at 1e-16. The one larger change is Spectral_PCC, at 4.3e-7 relative
(details below). The PR adds two new APIs and changes no existing signature.
Motivation: per FOV, the A549 ER eval on an A40 spends 7.2 s in MicroMS3IM scoring, 7.1 s in its
calibration, about 6 s in per-cell GLCM and 4.7 s in pixel metrics (dynacell
eval_timing.csvand py-spy).These follow the GPU segmentation in #57.
Changes
mean S,mean dS/dα) as one CuPyReductionKernelper evaluation, rows reduced in separate blockslo,hiadjacent floats); float32 elements never met the `MicroMS3IM.score_stack(gt, pred, degenerate=...): batched scoring with per-slice data ranges (ms_ssim's per-image path takes an(N,)range)cubic.feature.glcm_features_by_label: every region's GLCM in one direction-majorbincountper direction; properties over bounded label chunksbincounts (opt-in shared memory; out-of-range ids, host ids and sm < 60 take thebincountpath)compute_pixel_metrics: 403.6 → 299.5 ms per timepointThe fit and the scoring together, against
origin/mainon the same inputs: fit 45.9 → 1.48 s (31×), αidentical (3.2737876176834106).
Numerical equivalence
mainand to the production run's.glcm_featuresloop;glcm_featurestoo;reduce_poweron GPU forms |F|²as re² + im² rather than squaring |F|.
The fused kernels compute each per-element term with the array path's float operations in the same order,
compiled with
--fmad=false. Only the float64 summation order differs. Thebincountpath they replacealso added atomically, in an arbitrary order.
Tests
score_stack:scoreon CPU and GPU, across three batch sizes;glcm_features_by_label:glcm_featureswith each region's mask, in 2D/3D, on CPU/GPU, with shared and per-regionranges;
bincountpath for complex64 and complex128, plain, sectioned and power-only;Test plan
ruff check .,ruff format --check .,mypy --ignore-missing-imports cubic/pyteston an A40: 862 passed, 2 skipped (exit 0)/simplify, Copilot,/reviewand/code-reviewfindings addressed: one commit per finding, threads resolved🤖 Generated with Claude Code
https://claude.ai/code/session_011hEQ1TZJUbmhYKcC14RYHT
Summary by Sourcery
Accelerate GPU metric evaluation with fused reductions and batched scoring while adding batched MicroMS3IM and per-label GLCM APIs.
New Features:
Bug Fixes:
Enhancements:
Documentation:
Tests: