From 0f7a8bca9bc6e4ae76b517bca1dfa34a5c9f5166 Mon Sep 17 00:00:00 2001 From: Jonathan Wright Date: Sat, 14 Nov 2026 07:24:58 +0100 Subject: [PATCH 1/5] Document signed pixels: cut >= 0, negatives in the powder only The policy issue #9 asked for, as the code already does it (tested by test_dot's signed negative values and test_cut_range): the cut is never negative, so a sparse list holds no negative pixels; the powder sums every unmasked pixel as it is, so dark-subtracted CCD negatives integrate, and a Pilatus's gap and bad-pixel markers (-1, -2) belong in the mask. Co-Authored-By: Claude Opus 5.5 --- docs/index.md | 13 +++++++++++++ src/__init__.py | 8 +++++++- 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/docs/index.md b/docs/index.md index 3d3fc41..6456a5d 100644 --- a/docs/index.md +++ b/docs/index.md @@ -35,6 +35,19 @@ the free-threaded builds (the GPU decoder and the conveyor need Python 3). Decoding releases the GIL, and objects are not shared between threads: use one per thread. +## Pixel types and negative values + +Eiger data are unsigned (u8, u16, u32), and those are the types the fast paths are tuned for. +Signed pixels (i8, i16, i32, i64), as older Pilatus detectors and dark-subtracted CCD frames +have, are decoded the same way, with these rules: + +- The cut is never negative (a negative cut is an error), so the sparse pixel list holds only + pixels above a cut of zero or more: negative pixels never appear in it. +- The integration uses every unmasked pixel as it is, negative ones included: the powder is + the matrix times the masked frame. That is right for CCD data, where (intensity - dark) is + negative by noise. The negative markers a Pilatus writes for its gaps and bad pixels (-1, + -2) belong in the detector mask; unmasked, they would be summed into the powder. + ## More than sparse pixels A processor can limit how many pixels a frame returns. The pixels are chosen as diff --git a/src/__init__.py b/src/__init__.py index 9ca0dd0..89a9ddb 100644 --- a/src/__init__.py +++ b/src/__init__.py @@ -505,7 +505,8 @@ class chunk2sparse: cut is compared in the pixel type: for integer pixels a fractional cut is floored and a cut above the type's maximum selects nothing; a negative cut - is an error. + is an error. So with signed pixels (i8 .. i64) negative values are never + returned. c(chunk, cut) -> npixels, (values, indices) for one frame (the arrays are full size: slice [:npixels]); c.multi(chunks, cut) -> the same per @@ -731,6 +732,11 @@ class chunk2sparseCSC: howmany, thresholds: as chunk2sparse; they limit the pixels returned only: the powder always uses every unmasked pixel + Signed pixels: the cut is never negative, so negative pixels are not in the returned + list, but the powder sums every unmasked pixel as it is, negative ones included (CCD + frames with the dark subtracted). A Pilatus's negative markers for gaps and bad pixels + (-1, -2) must be masked to keep them out of the powder. + The array of the powder is the object's own and is overwritten by the next call. The object owns copies of the matrix and the mask, taken here (read-only): later changes to the arrays you passed do not reach it. The mask is folded into a raw matrix once, here: From c54814cf19da30e5f752ad27bd4c2a744e5f7417 Mon Sep 17 00:00:00 2001 From: Jonathan Wright Date: Sat, 14 Nov 2026 07:26:57 +0100 Subject: [PATCH 2/5] Signed pixels in the SIMD collects; i16 in the low-planes values i8, i16 and i32 pixels took the scalar collect whatever the pipeline said. Every SIMD collect (avx512cs, avx512f, avx2cs, sse2, vsx, neon) now takes a template argument S: the same kernel on the same bits with a signed compare (AVX-512 epi compares, no sign bias for SSE2/AVX2, the signed NEON/VSX compares, the scalar tails through bslz4_gt). collect_gt call them with S; collect_nz are the unsigned ones on the same bits (!= 0 is > 0 there). The unsigned instances are unchanged. i16 joins u8/u16 in the low-planes values: the cut is never negative, and a negative pixel sets the high planes, so its block takes kcb as a u16 block >= 256 does. i8 stays on kcb (its negative bytes would pass the fused collects' unsigned byte compare). Zen 4, sparsify fps (auto pipeline, before -> after): silicon i16 as recorded / zeroed <= 4 285 -> 769 / 339 -> 1507 silicon i8 as recorded / zeroed <= 4 323 -> 985 / 370 -> 2464 ceo2 i32 as recorded / zeroed <= 150 121 -> 277 / 216 -> 420 (i32 now within 7% of u32; i16 is held back by the high planes the -1 gap pixels set, which lz4 has to decode.) Tests: test_collect_tiers runs every value for i8/i16/i32 with -1 and minimum pixels; test_lowplanes_fused adds i16 frames with -1 pixels. x86 1552 passed; POWER9 1390 passed; aarch64 cross-build clean. Co-Authored-By: Claude Opus 5.5 --- src/_pipeline.py | 40 +++-- src/pipeline/pipeline.h | 2 +- src/pipeline/registry.c | 9 +- src/steps/collect.hpp | 302 +++++++++++++++++++++++++++++------ test/test_collect_tiers.py | 16 +- test/test_dot.py | 10 +- test/test_lowplanes_fused.py | 32 ++-- test/test_lz4zero.py | 4 +- 8 files changed, 324 insertions(+), 91 deletions(-) diff --git a/src/_pipeline.py b/src/_pipeline.py index 1bfe1e6..eacacc7 100644 --- a/src/_pipeline.py +++ b/src/_pipeline.py @@ -45,24 +45,24 @@ ), "untranspose": ( "auto", - ("lowplanes-vbmi", "u8/u16 pixels; AVX-512 VBMI, VBMI2 and GFNI"), - ("lowplanes-avx2", "u8/u16 pixels; AVX2"), - ("lowplanes-vsx", "u8/u16 pixels; POWER8 VSX"), - ("lowplanes-c", "u8/u16 pixels; portable C"), + ("lowplanes-vbmi", "u8/u16/i16 pixels; AVX-512 VBMI, VBMI2 and GFNI"), + ("lowplanes-avx2", "u8/u16/i16 pixels; AVX2"), + ("lowplanes-vsx", "u8/u16/i16 pixels; POWER8 VSX"), + ("lowplanes-c", "u8/u16/i16 pixels; portable C"), ("kcb", "kcb's bit transpose (SSE2/AVX2/AVX-512 dispatched inside)"), ("neon", "bitshuffle's NEON kernels; aarch64"), ("scalar", "bitshuffle's scalar kernels"), - ("lowplanes-avx512bw", "u8/u16 pixels; AVX-512 F/BW, for CPUs without VBMI and GFNI"), + ("lowplanes-avx512bw", "u8/u16/i16 pixels; AVX-512 F/BW, for CPUs without VBMI and GFNI"), ), "collect": ( "auto", - ("avx512cs", "u8/u16/u32 pixels; AVX-512 BW/VL and VBMI2"), - ("avx2cs", "u8/u16/u32 pixels; AVX2"), - ("vsx", "u8/u16/u32 pixels; POWER VSX"), - ("neon", "u8/u16/u32 pixels; aarch64 NEON"), - ("sse2", "u8/u16/u32 pixels; SSE2"), + ("avx512cs", "u8/u16/u32/i8/i16/i32 pixels; AVX-512 BW/VL and VBMI2"), + ("avx2cs", "u8/u16/u32/i8/i16/i32 pixels; AVX2"), + ("vsx", "u8/u16/u32/i8/i16/i32 pixels; POWER VSX"), + ("neon", "u8/u16/u32/i8/i16/i32 pixels; aarch64 NEON"), + ("sse2", "u8/u16/u32/i8/i16/i32 pixels; SSE2"), ("scalar", "any pixel type"), - ("avx512f", "u8/u16/u32 pixels; AVX-512 F/BW, for CPUs without VBMI2"), + ("avx512f", "u8/u16/u32/i8/i16/i32 pixels; AVX-512 F/BW, for CPUs without VBMI2"), ), "route": ( "auto", @@ -98,6 +98,12 @@ ("f", 4): "f32", ("f", 8): "f64"} +# the pixel dtypes of the SIMD collects (signed ones compare as signed) +SIMD_COLLECT_DTYPES = ("u8", "u16", "u32", "i8", "i16", "i32") +# the pixel dtypes of the low-planes untransposes (i16: a block with a negative pixel takes kcb) +LOWPLANES_DTYPES = ("u8", "u16", "i16") + + def dtype_suffix(dtype): dtype = np.dtype(dtype) try: @@ -177,14 +183,14 @@ def auto_step(step, suffix, codec, mask_all_valid): if step == "mask": return NAMES[step].index("none" if mask_all_valid else "pixel") if step == "untranspose": - if suffix in ("u8", "u16"): + if suffix in LOWPLANES_DTYPES: for n in ("lowplanes-vbmi", "lowplanes-avx512bw", "lowplanes-avx2", "lowplanes-vsx", "lowplanes-c"): if available(step, NAMES[step].index(n)) == 1: return NAMES[step].index(n) return NAMES[step].index("kcb") if step == "collect": - if suffix in ("u8", "u16", "u32"): + if suffix in SIMD_COLLECT_DTYPES: for n in ("avx512cs", "avx512f", "avx2cs", "vsx", "neon", "sse2"): if available(step, NAMES[step].index(n)) == 1: return NAMES[step].index(n) @@ -231,12 +237,12 @@ def resolve(pipeline, dtype, codec, mask, dot_auto=None, matrix=False): raise ValueError("pipeline['mask']='none': the mask has %d masked pixels (use 'pixel' or 'planes')" % int((np.asarray(mask) == 0).sum())) uname = NAMES["untranspose"][vals[UNTRANSPOSE]] - if uname.startswith("lowplanes") and suffix not in ("u8", "u16"): - raise ValueError("pipeline['untranspose']=%r is for u8/u16 pixels, not %s (use kcb)" + if uname.startswith("lowplanes") and suffix not in LOWPLANES_DTYPES: + raise ValueError("pipeline['untranspose']=%r is for u8/u16/i16 pixels, not %s (use kcb)" % (uname, np.dtype(dtype))) cname = NAMES["collect"][vals[COLLECT]] - if cname != "scalar" and suffix not in ("u8", "u16", "u32"): - raise ValueError("pipeline['collect']=%r is for u8/u16/u32 pixels, not %s (use scalar)" + if cname != "scalar" and suffix not in SIMD_COLLECT_DTYPES: + raise ValueError("pipeline['collect']=%r is for u8/u16/u32/i8/i16/i32 pixels, not %s (use scalar)" % (cname, np.dtype(dtype))) if matrix: if vals[DOT] == 0: diff --git a/src/pipeline/pipeline.h b/src/pipeline/pipeline.h index 71a58d7..abbcf70 100644 --- a/src/pipeline/pipeline.h +++ b/src/pipeline/pipeline.h @@ -75,7 +75,7 @@ enum { #define BSLZ4_IS_LOWPLANES(u) (((u) >= BSLZ4_UNTRANSPOSE_LOWPLANES_VBMI && (u) <= BSLZ4_UNTRANSPOSE_LOWPLANES_C) || \ (u) == BSLZ4_UNTRANSPOSE_LOWPLANES_AVX512BW) -/* collect: the SIMD values are u8/u16/u32 only */ +/* collect: the SIMD values take u8/u16/u32 and i8/i16/i32 (a signed compare) */ enum { BSLZ4_COLLECT_AVX512CS = 1, BSLZ4_COLLECT_AVX2CS = 2, diff --git a/src/pipeline/registry.c b/src/pipeline/registry.c index 614196c..eadc4e0 100644 --- a/src/pipeline/registry.c +++ b/src/pipeline/registry.c @@ -26,8 +26,9 @@ typedef char bslz4_value_slots_cover_every_step #define BSLZ4_ALL_DTYPES 0x3FFu #define BSLZ4_INT_DTYPES 0x0FFu /* u8 .. i64 */ #define BSLZ4_INT32_DTYPES 0x077u /* u8 u16 u32 i8 i16 i32: sums of 2**31 values fit int64 */ -#define BSLZ4_U8_U16_U32 ((1u << 0) | (1u << 1) | (1u << 2)) /* the SIMD collects */ -#define BSLZ4_U8_U16 ((1u << 0) | (1u << 1)) /* the low-planes values: u8, and u16 blocks < 256 */ +#define BSLZ4_SIMD_COLLECT_DTYPES 0x077u /* the SIMD collects: u8 u16 u32 i8 i16 i32 (signed: a signed compare) */ +#define BSLZ4_LOWPLANES_DTYPES ((1u << 0) | (1u << 1) | (1u << 5)) /* the low-planes values: u8, and + u16 / i16 blocks in 0..255 (a negative i16 sets the high planes) */ /* Each dot: the C entry (layout) it is reached through, the pixel dtypes it * takes, the powder element size (bytes; 0 = the pixel dtype, a packed @@ -177,9 +178,9 @@ int bslz4_resolve(int dtype, const uint16_t *pipeline, int layout, bslz4_pipe *B } if (!is_dot) p->step[BSLZ4_STEP_ROUTE] = p->step[BSLZ4_STEP_DOT] = 0; const uint32_t bit = 1u << dtype; - if (BSLZ4_IS_LOWPLANES(p->step[BSLZ4_STEP_UNTRANSPOSE]) && !(bit & BSLZ4_U8_U16)) + if (BSLZ4_IS_LOWPLANES(p->step[BSLZ4_STEP_UNTRANSPOSE]) && !(bit & BSLZ4_LOWPLANES_DTYPES)) return BSLZ4_ERR_DTYPE; - if (p->step[BSLZ4_STEP_COLLECT] != BSLZ4_COLLECT_SCALAR && !(bit & BSLZ4_U8_U16_U32)) + if (p->step[BSLZ4_STEP_COLLECT] != BSLZ4_COLLECT_SCALAR && !(bit & BSLZ4_SIMD_COLLECT_DTYPES)) return BSLZ4_ERR_DTYPE; if (is_dot) { const bslz4_dot_desc *d = &bslz4_dots[p->step[BSLZ4_STEP_DOT]]; diff --git a/src/steps/collect.hpp b/src/steps/collect.hpp index 11d109a..aa6c5d1 100644 --- a/src/steps/collect.hpp +++ b/src/steps/collect.hpp @@ -5,7 +5,8 @@ * sparse-route compaction (mask>0 & value!=0, called here with cut==0 -- * identical to !=0 for unsigned T). * - * Six SIMD collect values, u8/u16/u32 only: avx512cs (AVX-512 F+BW+VL+VBMI2, + * Six SIMD collect values, for u8/u16/u32 and i8/i16/i32 pixels (signed ones: the + * same kernels on the same bits with a signed compare, template argument S): avx512cs (AVX-512 F+BW+VL+VBMI2, * compress-store), avx512f (AVX-512 F+BW+VL without VBMI2: dword compress), * avx2cs (AVX2, table-driven compaction), sse2 (x86-64 ABI * baseline, always present), VSX (POWER8+, probe in @@ -18,7 +19,7 @@ * BSLZ4_HAVE_*_COLLECT guards below make this file inert there. * * BSLZ4_COLLECT_AVX512CS AVX-512 F+BW+VL+VBMI2: compress-store of the - * selected pixels and their addresses; u8/u16/u32. + * selected pixels and their addresses; u8/u16/u32 (and signed). * BSLZ4_COLLECT_AVX512F AVX-512 F+BW+VL (Skylake-SP, Cascade Lake: no * VBMI2): compare at the pixel width, compress 16 * pixels at a time as dwords; u8/u16/u32. @@ -27,8 +28,8 @@ * BSLZ4_COLLECT_VSX POWER8+ (needs -maltivec -mvsx); u8/u16/u32. * BSLZ4_COLLECT_NEON AArch64 ASIMD (ARMv8-A baseline); u8/u16/u32. * BSLZ4_COLLECT_SSE2 x86-64 baseline, always present; u8/u16/u32. - * BSLZ4_COLLECT_SCALAR the plain scalar loop, any dtype; every dtype - * other than u8/u16/u32 takes it whatever the value. + * BSLZ4_COLLECT_SCALAR the plain scalar loop, any dtype; u64, i64 and the + * floats take it whatever the value. * * The value actually used is pipeline[BSLZ4_STEP_COLLECT], validated by * bslz4_resolve (src/pipeline/registry.c). Availability is decided by @@ -50,6 +51,7 @@ #include "../pipeline/pipeline.h" #include +#include #include "c2py_amd64.h" #include "c2py_arm64.h" @@ -82,6 +84,43 @@ namespace bslz4 { +/* v > cut on the raw (unsigned) bits of a pixel, as a signed comparison when S: the SIMD + * collects take S for signed pixels (i8/i16/i32), working on the same bits as the unsigned + * ones (the values they store are the pixels' own). */ +template +inline bool bslz4_gt(U a, U b) { + typedef typename std::make_signed::type I; + return S ? (I) a > (I) b : a > b; +} + +#if BSLZ4_HAVE_VSX_COLLECT +/* the VSX compares, unsigned or (S) signed */ +template +inline __vector __bool char bslz4_vsx_gt(__vector unsigned char v, __vector unsigned char c) { + return S ? vec_cmpgt((__vector signed char) v, (__vector signed char) c) : vec_cmpgt(v, c); +} +template +inline int bslz4_vsx_any_gt(__vector unsigned char v, __vector unsigned char c) { + return S ? vec_any_gt((__vector signed char) v, (__vector signed char) c) : vec_any_gt(v, c); +} +template +inline __vector __bool short bslz4_vsx_gt(__vector unsigned short v, __vector unsigned short c) { + return S ? vec_cmpgt((__vector signed short) v, (__vector signed short) c) : vec_cmpgt(v, c); +} +template +inline int bslz4_vsx_any_gt(__vector unsigned short v, __vector unsigned short c) { + return S ? vec_any_gt((__vector signed short) v, (__vector signed short) c) : vec_any_gt(v, c); +} +template +inline __vector __bool int bslz4_vsx_gt(__vector unsigned int v, __vector unsigned int c) { + return S ? vec_cmpgt((__vector signed int) v, (__vector signed int) c) : vec_cmpgt(v, c); +} +template +inline int bslz4_vsx_any_gt(__vector unsigned int v, __vector unsigned int c) { + return S ? vec_any_gt((__vector signed int) v, (__vector signed int) c) : vec_any_gt(v, c); +} +#endif + /* A NULL mask means the mask step is none, or the mask was already applied * to the planes (every pixel is valid; src/pipeline/driver.c, and * src/steps/kernels.cpp passes `w->no_mask ? NULL : w->mask`), so the @@ -174,6 +213,7 @@ inline bool bslz4_neon_collect_capable() { * scratch (block_elems entries) both cover. The entries past the returned * count are scratch the caller never reads. */ #if BSLZ4_HAVE_VBMI_GFNI +template __attribute__((target("avx512f,avx512bw,avx512vl,avx512vbmi2,popcnt"))) inline int bslz4_collect_avx512cs_u16(const uint16_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, @@ -187,7 +227,7 @@ inline int bslz4_collect_avx512cs_u16(const uint16_t *BSLZ4_RESTRICT block, const __m512i sixteen = _mm512_set1_epi32(16); for (; j + 32 <= block_elems; j += 32) { const __m512i v = _mm512_loadu_si512((const void *) &block[j]); - __mmask32 m = _mm512_cmpgt_epu16_mask(v, vcut); + __mmask32 m = (S ? _mm512_cmpgt_epi16_mask(v, vcut) : _mm512_cmpgt_epu16_mask(v, vcut)); if (mask) m &= _mm256_cmpneq_epu8_mask(_mm256_loadu_si256((const __m256i *) &mask[j + i0]), _mm256_setzero_si256()); @@ -202,7 +242,7 @@ inline int bslz4_collect_avx512cs_u16(const uint16_t *BSLZ4_RESTRICT block, npx += nlo + __builtin_popcount((unsigned) mhi); } for (; j < block_elems; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -213,6 +253,7 @@ inline int bslz4_collect_avx512cs_u16(const uint16_t *BSLZ4_RESTRICT block, /* u32 version of collect avx512cs: vpcompressd is plain AVX-512F. Same store * bound as the u16 one, with 16 entries per group. */ +template __attribute__((target("avx512f,avx512bw,avx512vl,popcnt"))) inline int bslz4_collect_avx512cs_u32(const uint32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, @@ -225,7 +266,7 @@ inline int bslz4_collect_avx512cs_u32(const uint32_t *BSLZ4_RESTRICT block, const __m512i iota = _mm512_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15); for (; j + 16 <= block_elems; j += 16) { const __m512i v = _mm512_loadu_si512((const void *) &block[j]); - __mmask16 m = _mm512_cmpgt_epu32_mask(v, vcut); + __mmask16 m = (S ? _mm512_cmpgt_epi32_mask(v, vcut) : _mm512_cmpgt_epu32_mask(v, vcut)); if (mask) m &= _mm_cmpneq_epu8_mask(_mm_loadu_si128((const __m128i *) &mask[j + i0]), _mm_setzero_si128()); @@ -236,7 +277,7 @@ inline int bslz4_collect_avx512cs_u32(const uint32_t *BSLZ4_RESTRICT block, npx += __builtin_popcount((unsigned) m); } for (; j < block_elems; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -250,6 +291,7 @@ inline int bslz4_collect_avx512cs_u32(const uint32_t *BSLZ4_RESTRICT block, * reach at most 64 entries past the count: before group j at most j of this * block's pixels were selected, so a store ends before i0 + j + 64 <= * i0 + block_elems. */ +template __attribute__((target("avx512f,avx512bw,avx512vl,avx512vbmi2,popcnt"))) inline int bslz4_collect_avx512cs_u8(const uint8_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, @@ -263,7 +305,7 @@ inline int bslz4_collect_avx512cs_u8(const uint8_t *BSLZ4_RESTRICT block, const __m512i sixteen = _mm512_set1_epi32(16); for (; j + 64 <= block_elems; j += 64) { const __m512i v = _mm512_loadu_si512((const void *) &block[j]); - __mmask64 m = _mm512_cmpgt_epu8_mask(v, vcut); + __mmask64 m = (S ? _mm512_cmpgt_epi8_mask(v, vcut) : _mm512_cmpgt_epu8_mask(v, vcut)); if (mask) { const __m512i mv = _mm512_loadu_si512((const void *) &mask[j + i0]); m &= _mm512_test_epi8_mask(mv, mv); @@ -279,7 +321,7 @@ inline int bslz4_collect_avx512cs_u8(const uint8_t *BSLZ4_RESTRICT block, } } for (; j < block_elems; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -318,12 +360,12 @@ inline int bslz4_avx512f_emit16(__m512i v32, __mmask16 mq, size_t base, int npx, return npx + __builtin_popcount((unsigned) mq); } -template +template inline int bslz4_collect_tail(const T *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t j, size_t n, T cut, int npx, T *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { for (; j < n; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -332,6 +374,7 @@ inline int bslz4_collect_tail(const T *BSLZ4_RESTRICT block, const uint8_t *BSLZ return npx; } +template __attribute__((target("avx512f,avx512bw,popcnt"))) inline int bslz4_collect_avx512f_u8(const uint8_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, uint8_t cut, @@ -341,7 +384,7 @@ inline int bslz4_collect_avx512f_u8(const uint8_t *BSLZ4_RESTRICT block, const u const __m512i vcut = _mm512_set1_epi8((char) cut); for (; j + 64 <= n; j += 64) { const __m512i v = _mm512_loadu_si512((const void *) &block[j]); - __mmask64 m = _mm512_cmpgt_epu8_mask(v, vcut); + __mmask64 m = (S ? _mm512_cmpgt_epi8_mask(v, vcut) : _mm512_cmpgt_epu8_mask(v, vcut)); if (mask) { const __m512i mv = _mm512_loadu_si512((const void *) &mask[j + i0]); m &= _mm512_test_epi8_mask(mv, mv); @@ -365,9 +408,10 @@ inline int bslz4_collect_avx512f_u8(const uint8_t *BSLZ4_RESTRICT block, const u npx = bslz4_avx512f_emit16(_mm512_cvtepu8_epi32(_mm512_extracti32x4_epi32(v, 3)), (__mmask16) (m >> 48), j + i0 + 48, npx, out_vals, out_adr); } - return bslz4_collect_tail(block, mask, i0, j, n, cut, npx, out_vals, out_adr); + return bslz4_collect_tail(block, mask, i0, j, n, cut, npx, out_vals, out_adr); } +template __attribute__((target("avx512f,avx512bw,popcnt"))) inline int bslz4_collect_avx512f_u16(const uint16_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, uint16_t cut, @@ -377,7 +421,7 @@ inline int bslz4_collect_avx512f_u16(const uint16_t *BSLZ4_RESTRICT block, const const __m512i vcut = _mm512_set1_epi16((short) cut); for (; j + 32 <= n; j += 32) { const __m512i v = _mm512_loadu_si512((const void *) &block[j]); - __mmask32 m = _mm512_cmpgt_epu16_mask(v, vcut); + __mmask32 m = (S ? _mm512_cmpgt_epi16_mask(v, vcut) : _mm512_cmpgt_epu16_mask(v, vcut)); if (mask) { const __m256i mv = _mm256_loadu_si256((const __m256i *) (const void *) &mask[j + i0]); m &= ~(uint32_t) _mm256_movemask_epi8(_mm256_cmpeq_epi8(mv, _mm256_setzero_si256())); @@ -396,9 +440,10 @@ inline int bslz4_collect_avx512f_u16(const uint16_t *BSLZ4_RESTRICT block, const npx = bslz4_avx512f_emit16(_mm512_cvtepu16_epi32(_mm512_extracti64x4_epi64(v, 1)), (__mmask16) (m >> 16), j + i0 + 16, npx, out_vals, out_adr); } - return bslz4_collect_tail(block, mask, i0, j, n, cut, npx, out_vals, out_adr); + return bslz4_collect_tail(block, mask, i0, j, n, cut, npx, out_vals, out_adr); } +template __attribute__((target("avx512f,avx512bw,popcnt"))) inline int bslz4_collect_avx512f_u32(const uint32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, uint32_t cut, @@ -408,14 +453,14 @@ inline int bslz4_collect_avx512f_u32(const uint32_t *BSLZ4_RESTRICT block, const const __m512i vcut = _mm512_set1_epi32((int) cut); for (; j + 16 <= n; j += 16) { const __m512i v = _mm512_loadu_si512((const void *) &block[j]); - __mmask16 m = _mm512_cmpgt_epu32_mask(v, vcut); + __mmask16 m = (S ? _mm512_cmpgt_epi32_mask(v, vcut) : _mm512_cmpgt_epu32_mask(v, vcut)); if (mask) m &= (__mmask16) ~(uint32_t) _mm_movemask_epi8( _mm_cmpeq_epi8(_mm_loadu_si128((const __m128i *) (const void *) &mask[j + i0]), _mm_setzero_si128())); if (!m) continue; npx = bslz4_avx512f_emit16(v, m, j + i0, npx, out_vals, out_adr); } - return bslz4_collect_tail(block, mask, i0, j, n, cut, npx, out_vals, out_adr); + return bslz4_collect_tail(block, mask, i0, j, n, cut, npx, out_vals, out_adr); } #endif /* BSLZ4_HAVE_AVX512_COLLECT */ @@ -476,6 +521,7 @@ inline const uint64_t *bslz4_avx2cs_table() { return t.lane; } +template __attribute__((target("avx2,popcnt"))) inline int bslz4_collect_avx2cs_u16(const uint16_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, @@ -485,7 +531,7 @@ inline int bslz4_collect_avx2cs_u16(const uint16_t *BSLZ4_RESTRICT block, const uint64_t *lut = bslz4_avx2cs_table(); int npx = 0; size_t j = 0; - const __m256i sign_bias = _mm256_set1_epi16((short) 0x8000); + const __m256i sign_bias = _mm256_set1_epi16((short) (S ? 0 : 0x8000)); const __m256i vcut_b = _mm256_xor_si256(_mm256_set1_epi16((short) cut), sign_bias); const __m256i iota8 = _mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7); const __m128i one = _mm_set1_epi8(1); @@ -524,7 +570,7 @@ inline int bslz4_collect_avx2cs_u16(const uint16_t *BSLZ4_RESTRICT block, } } for (; j < block_elems; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -533,6 +579,7 @@ inline int bslz4_collect_avx2cs_u16(const uint16_t *BSLZ4_RESTRICT block, return npx; } +template __attribute__((target("avx2"))) inline int bslz4_collect_avx2cs_u32(const uint32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, @@ -541,7 +588,7 @@ inline int bslz4_collect_avx2cs_u32(const uint32_t *BSLZ4_RESTRICT block, uint32_t *BSLZ4_RESTRICT out_adr) { int npx = 0; size_t j = 0; - const __m256i sign_bias = _mm256_set1_epi32((int) 0x80000000u); + const __m256i sign_bias = _mm256_set1_epi32((int) (S ? 0 : 0x80000000u)); const __m256i vcut_b = _mm256_xor_si256(_mm256_set1_epi32((int) cut), sign_bias); for (; j + 8 <= block_elems; j += 8) { uint32_t mnz8 = mask @@ -563,7 +610,7 @@ inline int bslz4_collect_avx2cs_u32(const uint32_t *BSLZ4_RESTRICT block, } } for (; j < block_elems; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -575,6 +622,7 @@ inline int bslz4_collect_avx2cs_u32(const uint32_t *BSLZ4_RESTRICT block, /* u8 version of collect avx2cs: 32 pixels per step, the selection byte of * each 8 indexing the lane table as for u16 (vpermd the indices, pshufb the * bytes; stored full width, 8 entries past the count at most). */ +template __attribute__((target("avx2,popcnt"))) inline int bslz4_collect_avx2cs_u8(const uint8_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, @@ -584,7 +632,7 @@ inline int bslz4_collect_avx2cs_u8(const uint8_t *BSLZ4_RESTRICT block, const uint64_t *lut = bslz4_avx2cs_table(); int npx = 0; size_t j = 0; - const __m256i bias = _mm256_set1_epi8((char) 0x80); + const __m256i bias = _mm256_set1_epi8((char) (S ? 0 : 0x80)); const __m256i vcut_b = _mm256_xor_si256(_mm256_set1_epi8((char) cut), bias); const __m256i iota8 = _mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7); for (; j + 32 <= block_elems; j += 32) { @@ -618,7 +666,7 @@ inline int bslz4_collect_avx2cs_u8(const uint8_t *BSLZ4_RESTRICT block, } } for (; j < block_elems; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -642,6 +690,7 @@ inline int bslz4_collect_avx2cs_u8(const uint8_t *BSLZ4_RESTRICT block, * _mm_movemask_ps directly, 4 lanes, no packing. */ +template __attribute__((target("sse2"))) inline int bslz4_collect_sse2_u16(const uint16_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, @@ -650,7 +699,7 @@ inline int bslz4_collect_sse2_u16(const uint16_t *BSLZ4_RESTRICT block, uint32_t *BSLZ4_RESTRICT out_adr) { int npx = 0; size_t j = 0; - const __m128i sign_bias = _mm_set1_epi16((short) 0x8000); + const __m128i sign_bias = _mm_set1_epi16((short) (S ? 0 : 0x8000)); const __m128i vcut_b = _mm_xor_si128(_mm_set1_epi16((short) cut), sign_bias); for (; j + 8 <= block_elems; j += 8) { uint32_t mnz8 = mask @@ -672,7 +721,7 @@ inline int bslz4_collect_sse2_u16(const uint16_t *BSLZ4_RESTRICT block, } } for (; j < block_elems; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -681,6 +730,7 @@ inline int bslz4_collect_sse2_u16(const uint16_t *BSLZ4_RESTRICT block, return npx; } +template __attribute__((target("sse2"))) inline int bslz4_collect_sse2_u32(const uint32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, @@ -689,7 +739,7 @@ inline int bslz4_collect_sse2_u32(const uint32_t *BSLZ4_RESTRICT block, uint32_t *BSLZ4_RESTRICT out_adr) { int npx = 0; size_t j = 0; - const __m128i sign_bias = _mm_set1_epi32((int) 0x80000000u); + const __m128i sign_bias = _mm_set1_epi32((int) (S ? 0 : 0x80000000u)); const __m128i vcut_b = _mm_xor_si128(_mm_set1_epi32((int) cut), sign_bias); for (; j + 4 <= block_elems; j += 4) { uint32_t mnz4 = 0xFu; @@ -714,7 +764,7 @@ inline int bslz4_collect_sse2_u32(const uint32_t *BSLZ4_RESTRICT block, } } for (; j < block_elems; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -724,6 +774,7 @@ inline int bslz4_collect_sse2_u32(const uint32_t *BSLZ4_RESTRICT block, } /* u8 version of collect sse2: 16 pixels per compare, then the bit loop. */ +template __attribute__((target("sse2"))) inline int bslz4_collect_sse2_u8(const uint8_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, @@ -732,7 +783,7 @@ inline int bslz4_collect_sse2_u8(const uint8_t *BSLZ4_RESTRICT block, uint32_t *BSLZ4_RESTRICT out_adr) { int npx = 0; size_t j = 0; - const __m128i bias = _mm_set1_epi8((char) 0x80); + const __m128i bias = _mm_set1_epi8((char) (S ? 0 : 0x80)); const __m128i vcut_b = _mm_xor_si128(_mm_set1_epi8((char) cut), bias); for (; j + 16 <= block_elems; j += 16) { uint32_t mnz = mask @@ -750,7 +801,7 @@ inline int bslz4_collect_sse2_u8(const uint8_t *BSLZ4_RESTRICT block, } } for (; j < block_elems; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -807,6 +858,7 @@ static inline int bslz4_vsx_next_dense(int dense, int kept) { : (size_t) kept * BSLZ4_VSX_TO_DENSE > BSLZ4_VSX_WINDOW; } +template inline int bslz4_collect_vsx_u16(const uint16_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, uint16_t cut, @@ -822,7 +874,7 @@ inline int bslz4_collect_vsx_u16(const uint16_t *BSLZ4_RESTRICT block, const int start = npx; if (dense) { for (; j < end; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -832,8 +884,8 @@ inline int bslz4_collect_vsx_u16(const uint16_t *BSLZ4_RESTRICT block, for (; j < end; j += 8) { __vector unsigned short v; memcpy(&v, &block[j], sizeof(v)); - if (!vec_any_gt(v, vcut)) continue; - __vector __bool short gt = vec_cmpgt(v, vcut); + if (!bslz4_vsx_any_gt(v, vcut)) continue; + __vector __bool short gt = bslz4_vsx_gt(v, vcut); for (int lane = 0; lane < 8; lane++) { if (bslz4_mask_ok(mask, j + i0 + lane) && vec_extract((__vector unsigned short) gt, lane)) { out_vals[npx] = vec_extract(v, lane); @@ -846,7 +898,7 @@ inline int bslz4_collect_vsx_u16(const uint16_t *BSLZ4_RESTRICT block, dense = bslz4_vsx_next_dense(dense, npx - start); } for (; j < n; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -855,6 +907,7 @@ inline int bslz4_collect_vsx_u16(const uint16_t *BSLZ4_RESTRICT block, return npx; } +template inline int bslz4_collect_vsx_u32(const uint32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, uint32_t cut, @@ -870,7 +923,7 @@ inline int bslz4_collect_vsx_u32(const uint32_t *BSLZ4_RESTRICT block, const int start = npx; if (dense) { for (; j < end; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -880,8 +933,8 @@ inline int bslz4_collect_vsx_u32(const uint32_t *BSLZ4_RESTRICT block, for (; j < end; j += 4) { __vector unsigned int v; memcpy(&v, &block[j], sizeof(v)); - if (!vec_any_gt(v, vcut)) continue; - __vector __bool int gt = vec_cmpgt(v, vcut); + if (!bslz4_vsx_any_gt(v, vcut)) continue; + __vector __bool int gt = bslz4_vsx_gt(v, vcut); for (int lane = 0; lane < 4; lane++) { if (bslz4_mask_ok(mask, j + i0 + lane) && vec_extract((__vector unsigned int) gt, lane)) { out_vals[npx] = vec_extract(v, lane); @@ -894,7 +947,7 @@ inline int bslz4_collect_vsx_u32(const uint32_t *BSLZ4_RESTRICT block, dense = bslz4_vsx_next_dense(dense, npx - start); } for (; j < n; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -904,6 +957,7 @@ inline int bslz4_collect_vsx_u32(const uint32_t *BSLZ4_RESTRICT block, } /* u8 version of collect vsx: 16 lanes of bytes. */ +template inline int bslz4_collect_vsx_u8(const uint8_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, uint8_t cut, @@ -919,7 +973,7 @@ inline int bslz4_collect_vsx_u8(const uint8_t *BSLZ4_RESTRICT block, const int start = npx; if (dense) { for (; j < end; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -929,8 +983,8 @@ inline int bslz4_collect_vsx_u8(const uint8_t *BSLZ4_RESTRICT block, for (; j < end; j += 16) { __vector unsigned char v; memcpy(&v, &block[j], sizeof(v)); - if (!vec_any_gt(v, vcut)) continue; - __vector __bool char gt = vec_cmpgt(v, vcut); + if (!bslz4_vsx_any_gt(v, vcut)) continue; + __vector __bool char gt = bslz4_vsx_gt(v, vcut); for (int lane = 0; lane < 16; lane++) { if (bslz4_mask_ok(mask, j + i0 + lane) && vec_extract((__vector unsigned char) gt, lane)) { out_vals[npx] = vec_extract(v, lane); @@ -943,7 +997,7 @@ inline int bslz4_collect_vsx_u8(const uint8_t *BSLZ4_RESTRICT block, dense = bslz4_vsx_next_dense(dense, npx - start); } for (; j < n; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -973,6 +1027,7 @@ inline int bslz4_collect_vsx_u8(const uint8_t *BSLZ4_RESTRICT block, * of fighting that restriction with an unrolled constant-index ladder. */ +template inline int bslz4_collect_neon_u16(const uint16_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, uint16_t cut, @@ -983,7 +1038,7 @@ inline int bslz4_collect_neon_u16(const uint16_t *BSLZ4_RESTRICT block, uint16x8_t vcut = vdupq_n_u16(cut); for (; j + 8 <= n; j += 8) { uint16x8_t v = vld1q_u16(&block[j]); - uint16x8_t gt = vcgtq_u16(v, vcut); + uint16x8_t gt = (S ? vcgtq_s16(vreinterpretq_s16_u16(v), vreinterpretq_s16_u16(vcut)) : vcgtq_u16(v, vcut)); if (vmaxvq_u16(gt) == 0) continue; uint16_t vbuf[8], gbuf[8]; vst1q_u16(vbuf, v); @@ -997,7 +1052,7 @@ inline int bslz4_collect_neon_u16(const uint16_t *BSLZ4_RESTRICT block, } } for (; j < n; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -1006,6 +1061,7 @@ inline int bslz4_collect_neon_u16(const uint16_t *BSLZ4_RESTRICT block, return npx; } +template inline int bslz4_collect_neon_u32(const uint32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, uint32_t cut, @@ -1016,7 +1072,7 @@ inline int bslz4_collect_neon_u32(const uint32_t *BSLZ4_RESTRICT block, uint32x4_t vcut = vdupq_n_u32(cut); for (; j + 4 <= n; j += 4) { uint32x4_t v = vld1q_u32(&block[j]); - uint32x4_t gt = vcgtq_u32(v, vcut); + uint32x4_t gt = (S ? vcgtq_s32(vreinterpretq_s32_u32(v), vreinterpretq_s32_u32(vcut)) : vcgtq_u32(v, vcut)); if (vmaxvq_u32(gt) == 0) continue; uint32_t vbuf[4], gbuf[4]; vst1q_u32(vbuf, v); @@ -1030,7 +1086,7 @@ inline int bslz4_collect_neon_u32(const uint32_t *BSLZ4_RESTRICT block, } } for (; j < n; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -1040,6 +1096,7 @@ inline int bslz4_collect_neon_u32(const uint32_t *BSLZ4_RESTRICT block, } /* u8 version of collect neon: 16 lanes of bytes. */ +template inline int bslz4_collect_neon_u8(const uint8_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, uint8_t cut, @@ -1050,7 +1107,7 @@ inline int bslz4_collect_neon_u8(const uint8_t *BSLZ4_RESTRICT block, uint8x16_t vcut = vdupq_n_u8(cut); for (; j + 16 <= n; j += 16) { uint8x16_t v = vld1q_u8(&block[j]); - uint8x16_t gt = vcgtq_u8(v, vcut); + uint8x16_t gt = (S ? vcgtq_s8(vreinterpretq_s8_u8(v), vreinterpretq_s8_u8(vcut)) : vcgtq_u8(v, vcut)); if (vmaxvq_u8(gt) == 0) continue; uint8_t vbuf[16], gbuf[16]; vst1q_u8(vbuf, v); @@ -1064,7 +1121,7 @@ inline int bslz4_collect_neon_u8(const uint8_t *BSLZ4_RESTRICT block, } } for (; j < n; j++) { - if (bslz4_mask_ok(mask, j + i0) & (block[j] > cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -1342,6 +1399,155 @@ inline int bslz4_collect_nz(const uint8_t *BSLZ4_RESTRICT block, const return npx; } +/* signed pixels: the SIMD collects with a signed compare */ + +template<> +inline int bslz4_collect_gt(const int8_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, + size_t i0, size_t n, int8_t cut, int collect, + int8_t *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { + /* the unsigned kernels on the same bits, comparing as signed */ + const uint8_t *ub = (const uint8_t *) (const void *) block; + uint8_t *uo = (uint8_t *) (void *) out_vals; + const uint8_t uc = (uint8_t) cut; + (void) ub; (void) uo; (void) uc; +#if BSLZ4_HAVE_AVX512_COLLECT +#if BSLZ4_HAVE_VBMI_GFNI + if (collect == BSLZ4_COLLECT_AVX512CS) return bslz4_collect_avx512cs_u8(ub, mask, i0, n, uc, uo, out_adr); +#endif +#endif +#if BSLZ4_HAVE_AVX512_COLLECT + if (collect == BSLZ4_COLLECT_AVX512F) return bslz4_collect_avx512f_u8(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_AVX2_COLLECT + if (collect == BSLZ4_COLLECT_AVX2CS) return bslz4_collect_avx2cs_u8(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_SSE2_COLLECT + if (collect == BSLZ4_COLLECT_SSE2) return bslz4_collect_sse2_u8(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_VSX_COLLECT + if (collect == BSLZ4_COLLECT_VSX) return bslz4_collect_vsx_u8(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_NEON_COLLECT + if (collect == BSLZ4_COLLECT_NEON) return bslz4_collect_neon_u8(ub, mask, i0, n, uc, uo, out_adr); +#endif + int npx = 0; + for (size_t j = 0; j < n; j++) { + if (BSLZ4_UNLIKELY(bslz4_mask_ok(mask, j + i0) & (block[j] > cut))) { + out_vals[npx] = block[j]; + out_adr[npx] = (uint32_t) (j + i0); + npx++; + } + } + return npx; +} + +template<> +inline int bslz4_collect_nz(const int8_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, + size_t i0, size_t n, int collect, + int8_t *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { + /* != 0 is the same test on the unsigned bits (> 0 there): the unsigned kernels, cut 0 */ + return bslz4_collect_nz((const uint8_t *) (const void *) block, mask, i0, n, collect, + (uint8_t *) (void *) out_vals, out_adr); +} + +template<> +inline int bslz4_collect_gt(const int16_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, + size_t i0, size_t n, int16_t cut, int collect, + int16_t *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { + /* the unsigned kernels on the same bits, comparing as signed */ + const uint16_t *ub = (const uint16_t *) (const void *) block; + uint16_t *uo = (uint16_t *) (void *) out_vals; + const uint16_t uc = (uint16_t) cut; + (void) ub; (void) uo; (void) uc; +#if BSLZ4_HAVE_AVX512_COLLECT +#if BSLZ4_HAVE_VBMI_GFNI + if (collect == BSLZ4_COLLECT_AVX512CS) return bslz4_collect_avx512cs_u16(ub, mask, i0, n, uc, uo, out_adr); +#endif +#endif +#if BSLZ4_HAVE_AVX512_COLLECT + if (collect == BSLZ4_COLLECT_AVX512F) return bslz4_collect_avx512f_u16(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_AVX2_COLLECT + if (collect == BSLZ4_COLLECT_AVX2CS) return bslz4_collect_avx2cs_u16(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_SSE2_COLLECT + if (collect == BSLZ4_COLLECT_SSE2) return bslz4_collect_sse2_u16(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_VSX_COLLECT + if (collect == BSLZ4_COLLECT_VSX) return bslz4_collect_vsx_u16(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_NEON_COLLECT + if (collect == BSLZ4_COLLECT_NEON) return bslz4_collect_neon_u16(ub, mask, i0, n, uc, uo, out_adr); +#endif + int npx = 0; + for (size_t j = 0; j < n; j++) { + if (BSLZ4_UNLIKELY(bslz4_mask_ok(mask, j + i0) & (block[j] > cut))) { + out_vals[npx] = block[j]; + out_adr[npx] = (uint32_t) (j + i0); + npx++; + } + } + return npx; +} + +template<> +inline int bslz4_collect_nz(const int16_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, + size_t i0, size_t n, int collect, + int16_t *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { + /* != 0 is the same test on the unsigned bits (> 0 there): the unsigned kernels, cut 0 */ + return bslz4_collect_nz((const uint16_t *) (const void *) block, mask, i0, n, collect, + (uint16_t *) (void *) out_vals, out_adr); +} + +template<> +inline int bslz4_collect_gt(const int32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, + size_t i0, size_t n, int32_t cut, int collect, + int32_t *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { + /* the unsigned kernels on the same bits, comparing as signed */ + const uint32_t *ub = (const uint32_t *) (const void *) block; + uint32_t *uo = (uint32_t *) (void *) out_vals; + const uint32_t uc = (uint32_t) cut; + (void) ub; (void) uo; (void) uc; +#if BSLZ4_HAVE_AVX512_COLLECT +#if BSLZ4_HAVE_VBMI_GFNI + if (collect == BSLZ4_COLLECT_AVX512CS) return bslz4_collect_avx512cs_u32(ub, mask, i0, n, uc, uo, out_adr); +#endif +#endif +#if BSLZ4_HAVE_AVX512_COLLECT + if (collect == BSLZ4_COLLECT_AVX512F) return bslz4_collect_avx512f_u32(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_AVX2_COLLECT + if (collect == BSLZ4_COLLECT_AVX2CS) return bslz4_collect_avx2cs_u32(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_SSE2_COLLECT + if (collect == BSLZ4_COLLECT_SSE2) return bslz4_collect_sse2_u32(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_VSX_COLLECT + if (collect == BSLZ4_COLLECT_VSX) return bslz4_collect_vsx_u32(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_NEON_COLLECT + if (collect == BSLZ4_COLLECT_NEON) return bslz4_collect_neon_u32(ub, mask, i0, n, uc, uo, out_adr); +#endif + int npx = 0; + for (size_t j = 0; j < n; j++) { + if (BSLZ4_UNLIKELY(bslz4_mask_ok(mask, j + i0) & (block[j] > cut))) { + out_vals[npx] = block[j]; + out_adr[npx] = (uint32_t) (j + i0); + npx++; + } + } + return npx; +} + +template<> +inline int bslz4_collect_nz(const int32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, + size_t i0, size_t n, int collect, + int32_t *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { + /* != 0 is the same test on the unsigned bits (> 0 there): the unsigned kernels, cut 0 */ + return bslz4_collect_nz((const uint32_t *) (const void *) block, mask, i0, n, collect, + (uint32_t *) (void *) out_vals, out_adr); +} + #endif /* any tier available */ } /* namespace bslz4 */ diff --git a/test/test_collect_tiers.py b/test/test_collect_tiers.py index 59a5b38..9387b5b 100644 --- a/test/test_collect_tiers.py +++ b/test/test_collect_tiers.py @@ -1,6 +1,7 @@ """Every available collect value (avx512cs, avx2cs, vsx, neon, sse2, scalar) must give exactly numpy's selection: the values > cut of the unmasked -pixels, with their flat indices, in order. Run for u8, u16 and u32 over sparse, +pixels, with their flat indices, in order. Run for u8, u16, u32 and i8, i16, i32 (with +negative pixels, as a Pilatus writes gaps and bad pixels; cuts are never negative) over sparse, medium, dense and all-maximum data, several cuts, with a mask and with an all-ones mask (mask step 'none'), on a shape that leaves a tail; and the CSC powder must agree across collect values on both routes. @@ -45,7 +46,14 @@ def make(kind, dt, rng): a = rng.poisson(30, n) else: # "max" a = np.full(n, top, np.int64) - return np.minimum(a, top).astype(dt) + a = np.minimum(a, top) + if np.iinfo(dt).min < 0: # signed: some pixels -1 (gaps), some at the minimum + a[rng.random(n) < 0.02] = -1 + a[rng.random(n) < 0.002] = np.iinfo(dt).min + return a.astype(dt) + + +DTYPES = [np.uint8, np.uint16, np.uint32, np.int8, np.int16, np.int32] def chunks_of(frames, path): @@ -55,7 +63,7 @@ def chunks_of(frames, path): return [ds.id.read_direct_chunk((i, 0, 0))[1] for i in range(len(frames))], b.detect_codec(ds) -@pytest.mark.parametrize("dt", [np.uint8, np.uint16, np.uint32]) +@pytest.mark.parametrize("dt", DTYPES) @pytest.mark.parametrize("kind", ["sparse", "medium", "dense", "max"]) @pytest.mark.parametrize("masked", [True, False]) def test_every_tier_matches_numpy(tmp_path, dt, kind, masked): @@ -80,7 +88,7 @@ def test_every_tier_matches_numpy(tmp_path, dt, kind, masked): assert (vals[f, :n] == flat[want]).all(), (tier, kind, cut, f) -@pytest.mark.parametrize("dt", [np.uint8, np.uint16, np.uint32]) +@pytest.mark.parametrize("dt", DTYPES) @pytest.mark.parametrize("route", ["dense", "sparse"]) def test_powder_same_for_every_tier(tmp_path, dt, route): rng = np.random.default_rng(22) diff --git a/test/test_dot.py b/test/test_dot.py index 794df38..0ace7eb 100644 --- a/test/test_dot.py +++ b/test/test_dot.py @@ -210,8 +210,8 @@ def test_every_available_untranspose_matches(): for uname in usable: for dsname in chunks: dt, shp, chunklist = chunks[dsname] - if uname.startswith("lowplanes") and np.dtype(dt) not in (np.uint8, np.uint16): - with pytest.raises(ValueError, match="is for u8/u16 pixels"): + if uname.startswith("lowplanes") and np.dtype(dt) not in (np.uint8, np.uint16, np.int16): + with pytest.raises(ValueError, match="is for u8/u16/i16 pixels"): chunk2sparse(1 - ai.mask, dtype=np.dtype(dt), pipeline={"untranspose": uname}) continue bufs = [c for (filt, c) in chunklist] @@ -448,12 +448,12 @@ def P(c, **kw): def test_bad_collect_values_refused(): # the removed tiers ("avx512", "avx2") are unknown names; a SIMD collect - # needs u8/u16/u32 pixels; route/dot need a matrix + # needs u8/u16/u32/i8/i16/i32 pixels; route/dot need a matrix for bad in ("avx512", "avx2", len(_pipeline.NAMES["collect"])): with pytest.raises(ValueError, match="unknown"): chunk2sparse(1 - ai.mask, dtype=np.uint16, pipeline={"collect": bad}) - with pytest.raises(ValueError, match="u8/u16/u32"): - chunk2sparse(1 - ai.mask, dtype=np.int16, pipeline={"collect": "sse2"}) + with pytest.raises(ValueError, match="u8/u16/u32/i8/i16/i32"): + chunk2sparse(1 - ai.mask, dtype=np.int64, pipeline={"collect": "sse2"}) with pytest.raises(ValueError, match="no matrix"): chunk2sparse(1 - ai.mask, dtype=np.uint16, pipeline={"route": "dense"}) diff --git a/test/test_lowplanes_fused.py b/test/test_lowplanes_fused.py index 697b8bb..8e80364 100644 --- a/test/test_lowplanes_fused.py +++ b/test/test_lowplanes_fused.py @@ -1,5 +1,6 @@ -"""The low-planes untranspose's fused transpose + collect (u8 blocks, and u16 -blocks with the high byte-planes empty): the plain sparse output and the sparse-route CSC +"""The low-planes untranspose's fused transpose + collect (u8 blocks, and u16 or i16 +blocks with the high byte-planes empty; i16 frames also hold -1 pixels, as a Pilatus writes +gaps, whose blocks take kcb): the plain sparse output and the sparse-route CSC powder must equal numpy at every cut -- including cuts at and above 255, where nothing in such a block passes -- for every usable lowplanes-* value, with the mask applied per pixel and in the planes (and 'none' for an all-ones mask), @@ -48,6 +49,17 @@ def make(kind, rng): return a.astype(np.uint16) +def as_dtype(frames, dtype): + """The frames in dtype; i16: some pixels -1 (the blocks holding them fall back to kcb).""" + out = frames.astype(dtype) + if np.dtype(dtype).kind == "i": + out[:, ::17, 150] = -1 + return out + + +DTYPES = [np.uint16, np.uint8, np.int16] + + def chunks_of(frames, path, nelems=0): with h5py.File(str(path), "w") as f: ds = f.create_dataset("d", data=frames, chunks=(1,) + frames.shape[1:], @@ -66,12 +78,12 @@ def masks(): if n.startswith("lowplanes") and _pipeline.available("untranspose", i) == 1] -@pytest.mark.parametrize("dtype", [np.uint16, np.uint8]) +@pytest.mark.parametrize("dtype", DTYPES) @pytest.mark.parametrize("nelems", [0, 16384]) # default blocks; 32 kB of u16, 16 kB of u8 @pytest.mark.parametrize("kind", ["very sparse", "clustered", "dense <256", "max 255"]) @pytest.mark.parametrize("mname", ["ones", "holes"]) def test_sparse(tmp_path, kind, mname, nelems, dtype): - frames = make(kind, np.random.default_rng(17)).astype(dtype) + frames = as_dtype(make(kind, np.random.default_rng(17)), dtype) chunks, codec = chunks_of(frames, tmp_path / "s.h5", nelems) mask = masks()[mname] mps = ("pixel", "planes") + (("none",) if mname == "ones" else ()) @@ -90,11 +102,11 @@ def test_sparse(tmp_path, kind, mname, nelems, dtype): assert (vals[f, :n] == flat[want]).all(), tag -@pytest.mark.parametrize("dtype", [np.uint16, np.uint8]) +@pytest.mark.parametrize("dtype", DTYPES) @pytest.mark.parametrize("kind", ["very sparse", "clustered", "max 255"]) @pytest.mark.parametrize("mname", ["ones", "holes"]) def test_csc_sparse_route(tmp_path, kind, mname, dtype): - frames = make(kind, np.random.default_rng(18)).astype(dtype) + frames = as_dtype(make(kind, np.random.default_rng(18)), dtype) chunks, codec = chunks_of(frames, tmp_path / "c.h5") mask = masks()[mname] npix = SHAPE[0] * SHAPE[1] @@ -116,19 +128,19 @@ def test_csc_sparse_route(tmp_path, kind, mname, dtype): def test_u8_auto_takes_lowplanes(): - """The automatic untranspose for u8 is a low-planes value, as for u16.""" + """The automatic untranspose for u8 and i16 is a low-planes value, as for u16.""" mask = np.ones(SHAPE, np.uint8) - for dt in (np.uint8, np.uint16): + for dt in (np.uint8, np.uint16, np.int16): integ = b.chunk2sparse(mask, dtype=dt, codec=b.CODEC_LZ4) assert b.describe(integ.pipeline)["untranspose"] in LOWPLANES, np.dtype(dt).name -@pytest.mark.parametrize("dtype", [np.uint16, np.uint8]) +@pytest.mark.parametrize("dtype", DTYPES) @pytest.mark.parametrize("kind", ["clustered", "dense <256", "max 255"]) def test_csc_dense_route(tmp_path, kind, dtype): """The dense route builds the pixel block from the low planes (u16; u8 blocks take kcb): the powder and the sparse output of every low-planes value equal numpy's.""" - frames = make(kind, np.random.default_rng(19)).astype(dtype) + frames = as_dtype(make(kind, np.random.default_rng(19)), dtype) chunks, codec = chunks_of(frames, tmp_path / "d.h5") mask = masks()["holes"] npix = SHAPE[0] * SHAPE[1] diff --git a/test/test_lz4zero.py b/test/test_lz4zero.py index 18e26bf..21856e9 100644 --- a/test/test_lz4zero.py +++ b/test/test_lz4zero.py @@ -142,9 +142,9 @@ def test_powder_unchanged(tmp_path, dt, kind, route): def test_bad_values_refused(): - """lowplanes-* is u8 and u16 only; 'none' needs a mask with no zeros; the old + """lowplanes-* is u8, u16 and i16 only; 'none' needs a mask with no zeros; the old option names are not pipeline values.""" - with pytest.raises(ValueError, match="is for u8/u16 pixels"): + with pytest.raises(ValueError, match="is for u8/u16/i16 pixels"): b.chunk2sparse(MASK, dtype=np.uint32, pipeline={"untranspose": "lowplanes-c"}) with pytest.raises(ValueError, match="masked pixels"): b.chunk2sparse(MASK, dtype=np.uint16, pipeline={"mask": "none"}) From 6c2a84f9f94b68706d0dc2cb3fbdf82ab28efe7d Mon Sep 17 00:00:00 2001 From: Jonathan Wright Date: Sat, 14 Nov 2026 07:28:32 +0100 Subject: [PATCH 3/5] Docs: per-frame Pilatus markers are not supported Correction to 0f7a8bc: older Pilatus detectors change their -1/-2 markers from frame to frame, which a static mask cannot follow. Say so: only always-marked pixels can be masked; moving markers are summed as values. Co-Authored-By: Claude Opus 5.5 --- docs/index.md | 6 ++++-- src/__init__.py | 5 +++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/docs/index.md b/docs/index.md index 6456a5d..6079334 100644 --- a/docs/index.md +++ b/docs/index.md @@ -45,8 +45,10 @@ have, are decoded the same way, with these rules: pixels above a cut of zero or more: negative pixels never appear in it. - The integration uses every unmasked pixel as it is, negative ones included: the powder is the matrix times the masked frame. That is right for CCD data, where (intensity - dark) is - negative by noise. The negative markers a Pilatus writes for its gaps and bad pixels (-1, - -2) belong in the detector mask; unmasked, they would be summed into the powder. + negative by noise. +- Older Pilatus detectors mark gaps and bad pixels with negative values (-1, -2) that can + change from frame to frame. Pixels that are always marked can go in the detector mask, but + markers that move are not supported: they are summed into the powder as values. ## More than sparse pixels diff --git a/src/__init__.py b/src/__init__.py index 89a9ddb..f8befc8 100644 --- a/src/__init__.py +++ b/src/__init__.py @@ -734,8 +734,9 @@ class chunk2sparseCSC: Signed pixels: the cut is never negative, so negative pixels are not in the returned list, but the powder sums every unmasked pixel as it is, negative ones included (CCD - frames with the dark subtracted). A Pilatus's negative markers for gaps and bad pixels - (-1, -2) must be masked to keep them out of the powder. + frames with the dark subtracted). Negative markers that change from frame to frame (older + Pilatus detectors: -1, -2) are not supported: they are summed into the powder as values; + only the pixels always marked can go in the mask. The array of the powder is the object's own and is overwritten by the next call. The object owns copies of the matrix and the mask, taken here (read-only): later changes to From b87e01084ea57bc39c736ccddc64907cacb86a8d Mon Sep 17 00:00:00 2001 From: Jonathan Wright Date: Sat, 14 Nov 2026 08:17:37 +0100 Subject: [PATCH 4/5] negatives="list": negative pixels in the list, for per-frame markers Older Pilatus detectors mark gaps and bad pixels with -1/-2 that change from frame to frame, which a fixed mask cannot follow. chunk2sparse and chunk2sparseCSC (and plan.context) take negatives="drop" (default, as before: CCD frames with the dark subtracted) or "list": for signed integer pixels the pixel list also holds every negative pixel, whatever the cut, and the caller corrects the powder and its normalisation from it afterwards (docs/index.md, "Pixel types and negative values", now describing both cases). The powder is unchanged either way. It costs nothing in the kernels: with a cut >= 0, "v > cut or v < 0" is the unsigned compare on the same bits, so the signed collects run with S = false. The flag rides on the pipeline's collect value (BSLZ4_COLLECT_NEGATIVES, taken off by bslz4_resolve into bslz4_pipe.list_negatives and put on bslz4_work.collect by the driver); the dot paths' list filters and bslz4_cut_above_max honour it; collect_nz drops it before the unsigned dispatch. Refused for unsigned and float pixels (markers in an unsigned container are above any cut already) and with howmany. test_negatives.py: the list against numpy for i8/i16/i32/i64, every usable collect value, low-planes for i16, chunks and from_frames, a cut above the type's maximum; the powder and the documented correction on both routes for csc, csc-run, padded, bsb-csr, csc-nosplit, bsb-csr-nosplit and csc-nosplit-int; the refusals. x86 1561 passed; POWER9 1399 passed; aarch64 cross-build clean. Co-Authored-By: Claude Opus 5.5 --- docs/index.md | 35 +++++++--- src/__init__.py | 71 +++++++++++++++---- src/pipeline/driver.c | 2 +- src/pipeline/pipeline.h | 7 ++ src/pipeline/registry.c | 5 +- src/steps/collect.hpp | 81 +++++++++++++++------- src/steps/dot_csc.hpp | 9 +-- src/steps/dot_fixed.hpp | 2 +- src/steps/dot_padded.hpp | 2 +- test/test_negatives.py | 145 +++++++++++++++++++++++++++++++++++++++ 10 files changed, 299 insertions(+), 60 deletions(-) create mode 100644 test/test_negatives.py diff --git a/docs/index.md b/docs/index.md index 6079334..37a6e24 100644 --- a/docs/index.md +++ b/docs/index.md @@ -38,17 +38,30 @@ shared between threads: use one per thread. ## Pixel types and negative values Eiger data are unsigned (u8, u16, u32), and those are the types the fast paths are tuned for. -Signed pixels (i8, i16, i32, i64), as older Pilatus detectors and dark-subtracted CCD frames -have, are decoded the same way, with these rules: - -- The cut is never negative (a negative cut is an error), so the sparse pixel list holds only - pixels above a cut of zero or more: negative pixels never appear in it. -- The integration uses every unmasked pixel as it is, negative ones included: the powder is - the matrix times the masked frame. That is right for CCD data, where (intensity - dark) is - negative by noise. -- Older Pilatus detectors mark gaps and bad pixels with negative values (-1, -2) that can - change from frame to frame. Pixels that are always marked can go in the detector mask, but - markers that move are not supported: they are summed into the powder as values. +Signed pixels (i8, i16, i32, i64) are decoded the same way. The cut is never negative (a +negative cut is an error), and the integration always uses every unmasked pixel as it is, +negative ones included: the powder is the matrix times the masked frame. What differs is +whether negative pixels are returned in the pixel list, which a processor chooses when it is +made (`negatives`). + +Dark-subtracted CCD frames (the default, `negatives="drop"`). Negative values are noise and +belong in the integration as they are; the pixel list holds the pixels above the cut only. +Such frames are usually already decompressed, so the processor's `from_frames` call (the back +half of the pipeline, no decompression) is the one that applies. + +Detectors that mark gaps and bad pixels with negative values that change from frame to frame +(older Pilatus detectors write -1 and -2; `negatives="list"`). A fixed mask cannot follow +such markers, so the processor returns every negative pixel in the pixel list as well, at no +extra cost, and the caller corrects each frame afterwards from the list: with the marked +pixels' positions and values, subtract what they added to the powder (the matrix columns of +those pixels times their values), and subtract the same columns times the normalisation image +(solid angle, polarisation, flat field) from the normalisation sum. Both are sparse products +over a few columns, and the intensity is the corrected powder over the corrected +normalisation. The pixels always marked are better put in the mask. Markers stored in an +unsigned container (2**16-1, 2**32-1) need no option: they are above any cut, so they are in +the list already; the correction is the same, with the powder accumulated in float64 (or in +exact int64 sums for an integer matrix) so that subtracting such large values loses nothing +that matters. A limit on the pixels per frame (`howmany`) does not go with `negatives="list"`. ## More than sparse pixels diff --git a/src/__init__.py b/src/__init__.py index f8befc8..fb828ce 100644 --- a/src/__init__.py +++ b/src/__init__.py @@ -287,6 +287,34 @@ def _mask_auto_pending(pipeline, dtype, mask): and not bool((mask != 0).all())) +NEGATIVES_FLAG = 0x100 # BSLZ4_COLLECT_NEGATIVES (src/pipeline/pipeline.h) + + +def _negatives_flag(negatives, dtype, howmany): + """negatives="drop" (default) or "list": the flag OR'ed into the collect step for C.""" + if negatives in (None, "drop"): + return 0 + if negatives != "list": + raise ValueError('negatives must be "drop" or "list", not %r' % (negatives,)) + if np.dtype(dtype).kind != "i": + raise ValueError('negatives="list" is for signed integer pixels (i8, i16, i32, i64), not %s: ' + "in unsigned containers markers such as 2**32-1 are above any cut " + "already" % np.dtype(dtype)) + if howmany is not None: + raise ValueError('negatives="list" does not go with howmany (the thresholds would drop ' + "the negative pixels and the markers would use up the limit)") + return NEGATIVES_FLAG + + +def _with_negatives(pipe, flag): + """The pipeline array handed to C: the collect step with the negatives flag.""" + if not flag: + return pipe + p = pipe.copy() + p[_pipeline.COLLECT] |= flag + return p + + def _cut_array(cut, dtype): """The cut as a contiguous one-element array of the pixel dtype, which is how the C code takes it (it is compared in the pixel type). A Python number or a 0-d / one-element array. Integer @@ -502,11 +530,17 @@ class chunk2sparse: howmany of them in pixel order. thresholds = None: max(cut, 1) * 2**i for i in 0..15, or up to 16 sorted values. Only used when howmany is set. + negatives = "drop" (default) or "list", for signed integer pixels (i8 .. i64). "drop": the + list holds the pixels above the cut only, so negative pixels never appear in it (CCD + frames with the dark subtracted). "list": every negative pixel is in the list too, + whatever the cut -- for detectors that mark gaps and bad pixels with negative values + that change from frame to frame (older Pilatus: -1, -2), so the caller can correct + for them: see docs/index.md, "Pixel types and negative values". Not with howmany. cut is compared in the pixel type: for integer pixels a fractional cut is floored and a cut above the type's maximum selects nothing; a negative cut - is an error. So with signed pixels (i8 .. i64) negative values are never - returned. + is an error. So with signed pixels (i8 .. i64) negative values are only + returned with negatives="list". c(chunk, cut) -> npixels, (values, indices) for one frame (the arrays are full size: slice [:npixels]); c.multi(chunks, cut) -> the same per @@ -517,9 +551,12 @@ class chunk2sparse: """ def __init__( - self, mask, dtype=np.uint16, codec=None, pipeline=None, howmany=None, thresholds=None + self, mask, dtype=np.uint16, codec=None, pipeline=None, howmany=None, thresholds=None, + negatives="drop", ): self.howmany, self.thresholds = _toppix.check(howmany, thresholds) + self._negflag = _negatives_flag(negatives, dtype, howmany) + self.negatives = "list" if self._negflag else "drop" self.nfast = mask.shape[1] self.mask = mask.ravel() self.npix = mask.size @@ -585,7 +622,7 @@ def _run(self, pointers, lengths, cmp, cut, pixels=False): req.top, req.out, self._di, - pipe, + _with_negatives(pipe, self._negflag), ) if ret < 0: self._npx[:] = 0 @@ -729,14 +766,14 @@ class chunk2sparseCSC: pixels of at most 32 bits is auto "csc-nosplit-int": the sums are accumulated in int64 (no weights, no conversion to float) and returned as float64 like any float dot; a bin sum above 2**53 is rounded in that last cast - howmany, thresholds: as chunk2sparse; they limit the pixels returned - only: the powder always uses every unmasked pixel + howmany, thresholds, negatives: as chunk2sparse; they choose the pixels returned only: + the powder always uses every unmasked pixel as it is, negative ones included - Signed pixels: the cut is never negative, so negative pixels are not in the returned - list, but the powder sums every unmasked pixel as it is, negative ones included (CCD - frames with the dark subtracted). Negative markers that change from frame to frame (older - Pilatus detectors: -1, -2) are not supported: they are summed into the powder as values; - only the pixels always marked can go in the mask. + Signed pixels: the powder sums negative pixels as values either way (right for CCD frames + with the dark subtracted). For negative markers that change from frame to frame (older + Pilatus: -1, -2), negatives="list" returns them in the pixel list, and the caller removes + them from the powder and from the normalisation afterwards (docs/index.md, "Pixel types + and negative values"). The array of the powder is the object's own and is overwritten by the next call. The object owns copies of the matrix and the mask, taken here (read-only): later changes to @@ -758,8 +795,11 @@ def __init__( pipeline=None, howmany=None, thresholds=None, + negatives="drop", ): self.howmany, self.thresholds = _toppix.check(howmany, thresholds) + self._negflag = _negatives_flag(negatives, dtype, howmany) + self.negatives = "list" if self._negflag else "drop" # the matrix and mask are copies, so a later change to the caller's arrays changes nothing # here; a plan is used as it is (shared, its arrays read-only) if isinstance(csc, _matrix._Family): @@ -886,7 +926,7 @@ def _run(self, pointers, lengths, cmp, cut, pixels=False): req = _toppix.Request(self.howmany, self.thresholds, cut_a[0].item(), self.dtype, len(pointers)) head = (pointers, lengths, self.mask, self._values.ravel(), self._indices.ravel(), self._npx, cut_a, self._powder.ravel()) - tail = (self._di, pipe) + tail = (self._di, _with_negatives(pipe, self._negflag)) if lay.entry == "csc": ret = _ext.sparsify_and_dot( *(head + a + (self._workspace, self._cursors, req.top, req.out, lay.nout) + tail) @@ -1004,8 +1044,9 @@ def __init__(self, mask, csc, dtype=np.uint16, codec=None, pipeline=None): self.context() # refuses a pipeline or matrix that do not go together now, and makes # the first layout - def context(self, howmany=None, thresholds=None): + def context(self, howmany=None, thresholds=None, negatives="drop"): """A chunk2sparseCSC on this plan, for one thread: its own workspace, output arrays and - pipeline; the matrix and layouts are the plan's. howmany, thresholds: as chunk2sparse.""" + pipeline; the matrix and layouts are the plan's. howmany, thresholds, negatives: as + chunk2sparse.""" return chunk2sparseCSC(self.mask.reshape(self.shape), self, self.dtype, self._defaults[0], - self._defaults[1], howmany, thresholds) + self._defaults[1], howmany, thresholds, negatives) diff --git a/src/pipeline/driver.c b/src/pipeline/driver.c index 0da57f1..5ffe2da 100644 --- a/src/pipeline/driver.c +++ b/src/pipeline/driver.c @@ -260,7 +260,7 @@ int bslz4_driver_run(const int64_t *BSLZ4_RESTRICT compressed_ptrs, * block or tail. */ bslz4_work wbase = {0}; wbase.dtype = p->dtype; - wbase.collect = collect; + wbase.collect = collect | (p->list_negatives ? BSLZ4_COLLECT_NEGATIVES : 0); wbase.dot = dot; wbase.no_mask = emask == NULL; wbase.threshold = cut_i; diff --git a/src/pipeline/pipeline.h b/src/pipeline/pipeline.h index abbcf70..c1cb70b 100644 --- a/src/pipeline/pipeline.h +++ b/src/pipeline/pipeline.h @@ -76,6 +76,12 @@ enum { (u) == BSLZ4_UNTRANSPOSE_LOWPLANES_AVX512BW) /* collect: the SIMD values take u8/u16/u32 and i8/i16/i32 (a signed compare) */ + +/* negatives="list" (signed integer pixels): OR'ed into the pipeline's collect value by Python and + * into bslz4_work.collect by the driver -- the pixel list also holds every negative pixel (for + * markers such as a Pilatus's -1/-2 that change per frame; the caller corrects the powder and + * its normalisation from the list). Not a step value: bslz4_resolve takes it off. */ +#define BSLZ4_COLLECT_NEGATIVES 0x100 enum { BSLZ4_COLLECT_AVX512CS = 1, BSLZ4_COLLECT_AVX2CS = 2, @@ -158,6 +164,7 @@ typedef struct bslz4_pipe { uint16_t step[BSLZ4_NSTEPS]; int dtype; /* pixel dtype index 0..9: u8 u16 u32 u64 i8 i16 i32 i64 f32 f64 */ size_t elem_size; /* sizeof(pixel dtype) */ + int list_negatives; /* BSLZ4_COLLECT_NEGATIVES given (signed integer pixels only) */ } bslz4_pipe; /* top_pixels inside the block loop (ImageD11's rules, see src/_toppix.py). Each frame starts at diff --git a/src/pipeline/registry.c b/src/pipeline/registry.c index eadc4e0..2f04b87 100644 --- a/src/pipeline/registry.c +++ b/src/pipeline/registry.c @@ -165,8 +165,10 @@ static size_t bslz4_elem_size(int dtype) { int bslz4_resolve(int dtype, const uint16_t *pipeline, int layout, bslz4_pipe *BSLZ4_RESTRICT p) { if (dtype < 0 || dtype > 9) return BSLZ4_ERR_DTYPE; const int is_dot = layout >= 0; + const int negatives = (pipeline[BSLZ4_STEP_COLLECT] & BSLZ4_COLLECT_NEGATIVES) != 0; + if (negatives && !(dtype >= 4 && dtype <= 7)) return BSLZ4_ERR_DTYPE; /* i8 .. i64 only */ for (int s = 0; s < BSLZ4_NSTEPS; s++) { - const int v = pipeline[s]; + const int v = s == BSLZ4_STEP_COLLECT ? (pipeline[s] & ~BSLZ4_COLLECT_NEGATIVES) : pipeline[s]; if (!is_dot && (s == BSLZ4_STEP_ROUTE || s == BSLZ4_STEP_DOT)) { if (v != 0) return BSLZ4_ERR_BAD_PIPELINE; /* a plain sparsify has no dot */ continue; @@ -189,6 +191,7 @@ int bslz4_resolve(int dtype, const uint16_t *pipeline, int layout, bslz4_pipe *B } p->dtype = dtype; p->elem_size = bslz4_elem_size(dtype); + p->list_negatives = negatives; return 0; } diff --git a/src/steps/collect.hpp b/src/steps/collect.hpp index aa6c5d1..87e777f 100644 --- a/src/steps/collect.hpp +++ b/src/steps/collect.hpp @@ -93,6 +93,20 @@ inline bool bslz4_gt(U a, U b) { return S ? (I) a > (I) b : a > b; } +/* v < 0, false for unsigned and float-free of -Wtype-limits */ +template::value> struct bslz4_neg { + static inline bool of(T) { return false; } +}; +template struct bslz4_neg { + static inline bool of(T v) { return v < T(0); } +}; + +/* a pixel of the list: above the cut, or negative with negatives="list" */ +template +inline bool bslz4_keep(T v, T cut, int collect) { + return (v > cut) | (((collect & BSLZ4_COLLECT_NEGATIVES) != 0) & bslz4_neg::of(v)); +} + #if BSLZ4_HAVE_VSX_COLLECT /* the VSX compares, unsigned or (S) signed */ template @@ -1144,10 +1158,10 @@ template inline int bslz4_collect_gt(const T *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, T cut, int collect, T *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { - (void) collect; /* only the scalar tier is valid for non-SIMD dtypes */ + /* only the scalar tier is valid for non-SIMD dtypes; collect may carry BSLZ4_COLLECT_NEGATIVES */ int npx = 0; for (size_t j = 0; j < n; j++) { - if (BSLZ4_UNLIKELY(bslz4_mask_ok(mask, j + i0) & (block[j] > cut))) { + if (BSLZ4_UNLIKELY(bslz4_mask_ok(mask, j + i0) & bslz4_keep(block[j], cut, collect))) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -1410,29 +1424,33 @@ inline int bslz4_collect_gt(const int8_t *BSLZ4_RESTRICT block, const ui uint8_t *uo = (uint8_t *) (void *) out_vals; const uint8_t uc = (uint8_t) cut; (void) ub; (void) uo; (void) uc; + /* negatives="list": the unsigned compare on the same bits keeps v > cut and every v < 0 + * (the cut is never negative), so the same kernels serve with S = false */ + const bool neg = (collect & BSLZ4_COLLECT_NEGATIVES) != 0; + collect &= ~BSLZ4_COLLECT_NEGATIVES; #if BSLZ4_HAVE_AVX512_COLLECT #if BSLZ4_HAVE_VBMI_GFNI - if (collect == BSLZ4_COLLECT_AVX512CS) return bslz4_collect_avx512cs_u8(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_AVX512CS) return neg ? bslz4_collect_avx512cs_u8(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_avx512cs_u8(ub, mask, i0, n, uc, uo, out_adr); #endif #endif #if BSLZ4_HAVE_AVX512_COLLECT - if (collect == BSLZ4_COLLECT_AVX512F) return bslz4_collect_avx512f_u8(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_AVX512F) return neg ? bslz4_collect_avx512f_u8(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_avx512f_u8(ub, mask, i0, n, uc, uo, out_adr); #endif #if BSLZ4_HAVE_AVX2_COLLECT - if (collect == BSLZ4_COLLECT_AVX2CS) return bslz4_collect_avx2cs_u8(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_AVX2CS) return neg ? bslz4_collect_avx2cs_u8(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_avx2cs_u8(ub, mask, i0, n, uc, uo, out_adr); #endif #if BSLZ4_HAVE_SSE2_COLLECT - if (collect == BSLZ4_COLLECT_SSE2) return bslz4_collect_sse2_u8(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_SSE2) return neg ? bslz4_collect_sse2_u8(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_sse2_u8(ub, mask, i0, n, uc, uo, out_adr); #endif #if BSLZ4_HAVE_VSX_COLLECT - if (collect == BSLZ4_COLLECT_VSX) return bslz4_collect_vsx_u8(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_VSX) return neg ? bslz4_collect_vsx_u8(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_vsx_u8(ub, mask, i0, n, uc, uo, out_adr); #endif #if BSLZ4_HAVE_NEON_COLLECT - if (collect == BSLZ4_COLLECT_NEON) return bslz4_collect_neon_u8(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_NEON) return neg ? bslz4_collect_neon_u8(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_neon_u8(ub, mask, i0, n, uc, uo, out_adr); #endif int npx = 0; for (size_t j = 0; j < n; j++) { - if (BSLZ4_UNLIKELY(bslz4_mask_ok(mask, j + i0) & (block[j] > cut))) { + if (BSLZ4_UNLIKELY(bslz4_mask_ok(mask, j + i0) & ((block[j] > cut) | (neg & (block[j] < 0))))) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -1446,7 +1464,8 @@ inline int bslz4_collect_nz(const int8_t *BSLZ4_RESTRICT block, const ui size_t i0, size_t n, int collect, int8_t *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { /* != 0 is the same test on the unsigned bits (> 0 there): the unsigned kernels, cut 0 */ - return bslz4_collect_nz((const uint8_t *) (const void *) block, mask, i0, n, collect, + return bslz4_collect_nz((const uint8_t *) (const void *) block, mask, i0, n, + collect & ~BSLZ4_COLLECT_NEGATIVES, (uint8_t *) (void *) out_vals, out_adr); } @@ -1459,29 +1478,33 @@ inline int bslz4_collect_gt(const int16_t *BSLZ4_RESTRICT block, const uint16_t *uo = (uint16_t *) (void *) out_vals; const uint16_t uc = (uint16_t) cut; (void) ub; (void) uo; (void) uc; + /* negatives="list": the unsigned compare on the same bits keeps v > cut and every v < 0 + * (the cut is never negative), so the same kernels serve with S = false */ + const bool neg = (collect & BSLZ4_COLLECT_NEGATIVES) != 0; + collect &= ~BSLZ4_COLLECT_NEGATIVES; #if BSLZ4_HAVE_AVX512_COLLECT #if BSLZ4_HAVE_VBMI_GFNI - if (collect == BSLZ4_COLLECT_AVX512CS) return bslz4_collect_avx512cs_u16(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_AVX512CS) return neg ? bslz4_collect_avx512cs_u16(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_avx512cs_u16(ub, mask, i0, n, uc, uo, out_adr); #endif #endif #if BSLZ4_HAVE_AVX512_COLLECT - if (collect == BSLZ4_COLLECT_AVX512F) return bslz4_collect_avx512f_u16(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_AVX512F) return neg ? bslz4_collect_avx512f_u16(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_avx512f_u16(ub, mask, i0, n, uc, uo, out_adr); #endif #if BSLZ4_HAVE_AVX2_COLLECT - if (collect == BSLZ4_COLLECT_AVX2CS) return bslz4_collect_avx2cs_u16(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_AVX2CS) return neg ? bslz4_collect_avx2cs_u16(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_avx2cs_u16(ub, mask, i0, n, uc, uo, out_adr); #endif #if BSLZ4_HAVE_SSE2_COLLECT - if (collect == BSLZ4_COLLECT_SSE2) return bslz4_collect_sse2_u16(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_SSE2) return neg ? bslz4_collect_sse2_u16(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_sse2_u16(ub, mask, i0, n, uc, uo, out_adr); #endif #if BSLZ4_HAVE_VSX_COLLECT - if (collect == BSLZ4_COLLECT_VSX) return bslz4_collect_vsx_u16(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_VSX) return neg ? bslz4_collect_vsx_u16(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_vsx_u16(ub, mask, i0, n, uc, uo, out_adr); #endif #if BSLZ4_HAVE_NEON_COLLECT - if (collect == BSLZ4_COLLECT_NEON) return bslz4_collect_neon_u16(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_NEON) return neg ? bslz4_collect_neon_u16(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_neon_u16(ub, mask, i0, n, uc, uo, out_adr); #endif int npx = 0; for (size_t j = 0; j < n; j++) { - if (BSLZ4_UNLIKELY(bslz4_mask_ok(mask, j + i0) & (block[j] > cut))) { + if (BSLZ4_UNLIKELY(bslz4_mask_ok(mask, j + i0) & ((block[j] > cut) | (neg & (block[j] < 0))))) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -1495,7 +1518,8 @@ inline int bslz4_collect_nz(const int16_t *BSLZ4_RESTRICT block, const size_t i0, size_t n, int collect, int16_t *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { /* != 0 is the same test on the unsigned bits (> 0 there): the unsigned kernels, cut 0 */ - return bslz4_collect_nz((const uint16_t *) (const void *) block, mask, i0, n, collect, + return bslz4_collect_nz((const uint16_t *) (const void *) block, mask, i0, n, + collect & ~BSLZ4_COLLECT_NEGATIVES, (uint16_t *) (void *) out_vals, out_adr); } @@ -1508,29 +1532,33 @@ inline int bslz4_collect_gt(const int32_t *BSLZ4_RESTRICT block, const uint32_t *uo = (uint32_t *) (void *) out_vals; const uint32_t uc = (uint32_t) cut; (void) ub; (void) uo; (void) uc; + /* negatives="list": the unsigned compare on the same bits keeps v > cut and every v < 0 + * (the cut is never negative), so the same kernels serve with S = false */ + const bool neg = (collect & BSLZ4_COLLECT_NEGATIVES) != 0; + collect &= ~BSLZ4_COLLECT_NEGATIVES; #if BSLZ4_HAVE_AVX512_COLLECT #if BSLZ4_HAVE_VBMI_GFNI - if (collect == BSLZ4_COLLECT_AVX512CS) return bslz4_collect_avx512cs_u32(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_AVX512CS) return neg ? bslz4_collect_avx512cs_u32(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_avx512cs_u32(ub, mask, i0, n, uc, uo, out_adr); #endif #endif #if BSLZ4_HAVE_AVX512_COLLECT - if (collect == BSLZ4_COLLECT_AVX512F) return bslz4_collect_avx512f_u32(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_AVX512F) return neg ? bslz4_collect_avx512f_u32(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_avx512f_u32(ub, mask, i0, n, uc, uo, out_adr); #endif #if BSLZ4_HAVE_AVX2_COLLECT - if (collect == BSLZ4_COLLECT_AVX2CS) return bslz4_collect_avx2cs_u32(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_AVX2CS) return neg ? bslz4_collect_avx2cs_u32(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_avx2cs_u32(ub, mask, i0, n, uc, uo, out_adr); #endif #if BSLZ4_HAVE_SSE2_COLLECT - if (collect == BSLZ4_COLLECT_SSE2) return bslz4_collect_sse2_u32(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_SSE2) return neg ? bslz4_collect_sse2_u32(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_sse2_u32(ub, mask, i0, n, uc, uo, out_adr); #endif #if BSLZ4_HAVE_VSX_COLLECT - if (collect == BSLZ4_COLLECT_VSX) return bslz4_collect_vsx_u32(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_VSX) return neg ? bslz4_collect_vsx_u32(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_vsx_u32(ub, mask, i0, n, uc, uo, out_adr); #endif #if BSLZ4_HAVE_NEON_COLLECT - if (collect == BSLZ4_COLLECT_NEON) return bslz4_collect_neon_u32(ub, mask, i0, n, uc, uo, out_adr); + if (collect == BSLZ4_COLLECT_NEON) return neg ? bslz4_collect_neon_u32(ub, mask, i0, n, uc, uo, out_adr) : bslz4_collect_neon_u32(ub, mask, i0, n, uc, uo, out_adr); #endif int npx = 0; for (size_t j = 0; j < n; j++) { - if (BSLZ4_UNLIKELY(bslz4_mask_ok(mask, j + i0) & (block[j] > cut))) { + if (BSLZ4_UNLIKELY(bslz4_mask_ok(mask, j + i0) & ((block[j] > cut) | (neg & (block[j] < 0))))) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -1544,7 +1572,8 @@ inline int bslz4_collect_nz(const int32_t *BSLZ4_RESTRICT block, const size_t i0, size_t n, int collect, int32_t *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { /* != 0 is the same test on the unsigned bits (> 0 there): the unsigned kernels, cut 0 */ - return bslz4_collect_nz((const uint32_t *) (const void *) block, mask, i0, n, collect, + return bslz4_collect_nz((const uint32_t *) (const void *) block, mask, i0, n, + collect & ~BSLZ4_COLLECT_NEGATIVES, (uint32_t *) (void *) out_vals, out_adr); } diff --git a/src/steps/dot_csc.hpp b/src/steps/dot_csc.hpp index 91ad5f3..2ac39cb 100644 --- a/src/steps/dot_csc.hpp +++ b/src/steps/dot_csc.hpp @@ -26,7 +26,8 @@ using namespace bslz4; * "no sparse output" arrives as exactly the maximum. uint64_t and the * floating types are never skipped. */ template inline bool bslz4_cut_above_max(const bslz4_work *w) { - return w->threshold >= (int64_t) std::numeric_limits::max(); + return w->threshold >= (int64_t) std::numeric_limits::max() && + !(w->collect & BSLZ4_COLLECT_NEGATIVES); /* negatives are still listed */ } /* The cut as the pixel type: integers from the int64 (u64: its bit pattern), floats from cutf, so * a fractional cut is exact. */ @@ -184,7 +185,7 @@ int sparse_dot_core(const bslz4_work *BSLZ4_RESTRICT w, const bslz4_mat_csc *BSL const int skip = bslz4_cut_above_max(w); for (int kk = 0; kk < nz; kk++) { T val = tv[kk]; - if (!skip && BSLZ4_UNLIKELY(val > cut)) { + if (!skip && BSLZ4_UNLIKELY(bslz4_keep(val, cut, w->collect))) { ov[npx] = val; out_adr[npx] = tidx[kk]; npx++; @@ -269,7 +270,7 @@ int run_sparse(const bslz4_work *BSLZ4_RESTRICT w) { const int skip = bslz4_cut_above_max(w); for (int kk = 0; kk < nz; kk++) { T val = tv[kk]; - if (!skip && BSLZ4_UNLIKELY(val > cut)) { + if (!skip && BSLZ4_UNLIKELY(bslz4_keep(val, cut, w->collect))) { ov[npx] = val; out_adr[npx] = tidx[kk]; npx++; @@ -318,7 +319,7 @@ int nosplit_sparse_core(const bslz4_work *BSLZ4_RESTRICT w, const uint32_t *BSLZ const int skip = bslz4_cut_above_max(w); for (int kk = 0; kk < nz; kk++) { T val = tv[kk]; - if (!skip && BSLZ4_UNLIKELY(val > cut)) { + if (!skip && BSLZ4_UNLIKELY(bslz4_keep(val, cut, w->collect))) { ov[npx] = val; out_adr[npx] = tidx[kk]; npx++; diff --git a/src/steps/dot_fixed.hpp b/src/steps/dot_fixed.hpp index 0d07ea2..2dd3af8 100644 --- a/src/steps/dot_fixed.hpp +++ b/src/steps/dot_fixed.hpp @@ -181,7 +181,7 @@ static inline int sparse(const bslz4_work *BSLZ4_RESTRICT w) { const int skip = bslz4_cut_above_max(w); for (int kk = 0; kk < nz; kk++) { T val = tv[kk]; - if (!skip && BSLZ4_UNLIKELY(val > cut)) { + if (!skip && BSLZ4_UNLIKELY(bslz4_keep(val, cut, w->collect))) { ov[npx] = val; out_adr[npx] = tidx[kk]; npx++; diff --git a/src/steps/dot_padded.hpp b/src/steps/dot_padded.hpp index 756f790..ee9f762 100644 --- a/src/steps/dot_padded.hpp +++ b/src/steps/dot_padded.hpp @@ -193,7 +193,7 @@ int padded_sparse(const bslz4_work *BSLZ4_RESTRICT w) { const int skip = bslz4_cut_above_max(w); for (int kk = 0; kk < nz; kk++) { T val = tv[kk]; - if (!skip && BSLZ4_UNLIKELY(val > cut)) { + if (!skip && BSLZ4_UNLIKELY(bslz4_keep(val, cut, w->collect))) { ov[npx] = val; out_adr[npx] = tidx[kk]; npx++; diff --git a/test/test_negatives.py b/test/test_negatives.py new file mode 100644 index 0000000..958a410 --- /dev/null +++ b/test/test_negatives.py @@ -0,0 +1,145 @@ +"""negatives="list": for signed integer pixels the pixel list also holds every negative pixel, +whatever the cut (markers such as an older Pilatus's -1/-2 that change from frame to frame), +while the powder is unchanged (every unmasked pixel as it is). Checked against numpy for i8, +i16, i32 and i64, every usable collect value, the low-planes untranspose (i16), chunks and +uncompressed frames, both dot routes and several dots; the correction the docs describe +(docs/index.md, "Pixel types and negative values") recovers the powder of the unmarked pixels; +the option is refused where it does not apply. +""" + +import os +import sys + +_path = os.environ.get("BSLZ4_TO_SPARSE_PATH") +if _path: + sys.path.insert(0, _path) + +import numpy as np +import pytest + +h5py = pytest.importorskip("h5py") +hdf5plugin = pytest.importorskip("hdf5plugin") +sp = pytest.importorskip("scipy.sparse") + +import bslz4_to_sparse as b +from bslz4_to_sparse import _pipeline + +SHAPE = (93, 217) # full blocks plus a tail +DTYPES = [np.int8, np.int16, np.int32, np.int64] +COLLECTS = [n for i, n in enumerate(_pipeline.NAMES["collect"]) + if i and _pipeline.available("collect", i) == 1] +LOWPLANES = [n for i, n in enumerate(_pipeline.NAMES["untranspose"]) + if n.startswith("lowplanes") and _pipeline.available("untranspose", i) == 1] + + +def frames_with_markers(dt, rng, n=3): + """Poisson counts with -1/-2 markers at positions that change per frame.""" + f = np.minimum(rng.poisson(0.3, (n,) + SHAPE), 100).astype(np.int64) + for i in range(n): + d = rng.random(SHAPE) < 0.01 + f[i][d] = np.where(rng.random(int(d.sum())) < 0.5, -1, -2) + f[:, 5, 10:80] = 50 # a run above any small cut + return f.astype(dt) + + +def chunks_of(frames, path): + with h5py.File(str(path), "w") as h: + ds = h.create_dataset("d", data=frames, chunks=(1,) + frames.shape[1:], + **hdf5plugin.Bitshuffle(nelems=0, cname="lz4")) + return [ds.id.read_direct_chunk((i, 0, 0))[1] for i in range(len(frames))] + + +def mask_of(rng): + m = (rng.random(SHAPE) > 0.05).astype(np.uint8) + m[:, 100:103] = 0 + return m + + +def check_list(frames, mask, cut, npx, vals, adr, tag): + for f in range(len(frames)): + flat = frames[f].ravel() + want = np.flatnonzero(((flat > cut) | (flat < 0)) & (mask.ravel() > 0)) + n = int(npx[f]) + assert n == want.size, tag + (f, n, want.size) + assert (adr[f, :n] == want).all(), tag + (f,) + assert (vals[f, :n] == flat[want]).all(), tag + (f,) + + +@pytest.mark.parametrize("dt", DTYPES) +def test_list_holds_negatives(tmp_path, dt): + rng = np.random.default_rng(31) + frames = frames_with_markers(dt, rng) + chunks = chunks_of(frames, tmp_path / "n.h5") + mask = mask_of(rng) + big = int(np.iinfo(dt).max) + 10 + unt = ["kcb"] + (LOWPLANES if dt == np.int16 else []) + for collect in (COLLECTS if np.dtype(dt).itemsize < 8 else ["scalar"]): # i64: scalar only + for u in unt: + o = b.chunk2sparse(mask, dtype=dt, codec=b.CODEC_LZ4, negatives="list", + pipeline={"collect": collect, "untranspose": u}) + for cut in (0, 3, big): # above the maximum: the negatives only + npx, (v, a) = o.multi(chunks, cut) + check_list(frames, mask, min(cut, int(np.iinfo(dt).max)), npx, v, a, + ("chunks", np.dtype(dt).name, collect, u, cut)) + npx, (v, a) = o.from_frames(frames, cut) + check_list(frames, mask, min(cut, int(np.iinfo(dt).max)), npx, v, a, + ("frames", np.dtype(dt).name, collect, u, cut)) + # the default leaves them out, as before + d = b.chunk2sparse(mask, dtype=dt, codec=b.CODEC_LZ4, + pipeline={"collect": collect, "untranspose": u}) + npx, (v, a) = d.multi(chunks, 0) + assert all((v[f, :npx[f]] > 0).all() for f in range(len(frames))) + + +@pytest.mark.parametrize("dt", [np.int16, np.int32]) +@pytest.mark.parametrize("route", ["dense", "sparse"]) +def test_powder_and_correction(tmp_path, dt, route): + rng = np.random.default_rng(32) + frames = frames_with_markers(dt, rng) + chunks = chunks_of(frames, tmp_path / "p.h5") + mask = mask_of(rng) + npix = mask.size + M = sp.random(150, npix, density=0.02, format="csc", random_state=4, dtype=np.float32) + b0 = rng.integers(0, 149, npix) # each pixel in two consecutive bins + runs = sp.csc_matrix((np.tile([0.6, 0.4], npix).astype(np.float32), + np.stack([b0, b0 + 1], 1).ravel(), np.arange(0, 2 * npix + 1, 2)), + shape=(150, npix)) + nosplit = sp.csc_matrix((np.ones(npix, np.float32), rng.integers(0, 40, npix), + np.arange(npix + 1)), shape=(40, npix)) + norm = (0.5 + rng.random(npix)).astype(np.float32) # solid angle, polarisation, ... + cases = [(M, "csc")] + [(runs, d) for d in ("csc-run", "padded", "bsb-csr")] + \ + [(nosplit, d) for d in ("csc-nosplit", "bsb-csr-nosplit", "csc-nosplit-int")] + for mat, dot in cases: + for collect in COLLECTS: + o = b.chunk2sparseCSC(mask, mat, dtype=dt, codec=b.CODEC_LZ4, negatives="list", + pipeline={"route": route, "dot": dot, "collect": collect}) + npx, (v, a), p = o.multi(chunks, 2) + check_list(frames, mask, 2, npx, v, a, (np.dtype(dt).name, route, dot, collect)) + Md = mat.tocsc().astype(np.float64) + for f in range(len(frames)): + flat = frames[f].ravel().astype(np.float64) + m = mask.ravel() > 0 + assert np.allclose(p[f], Md @ (flat * m), rtol=1e-9, atol=1e-9), (dot, route, f) + # the caller's correction from the list: the marked pixels out of the powder and + # out of the normalisation + n = int(npx[f]) + vals, idx = v[f, :n], a[f, :n].astype(np.int64) + marked = vals < 0 + fixed = p[f] - Md[:, idx[marked]] @ vals[marked].astype(np.float64) + keep = m & (flat >= 0) + assert np.allclose(fixed, Md @ (flat * keep), rtol=1e-9, atol=1e-9), (dot, route, f) + nfix = Md @ (norm * m) - Md[:, idx[marked]] @ norm[idx[marked]] + assert np.allclose(nfix, Md @ (norm * keep), rtol=1e-6, atol=1e-6), (dot, route, f) + + +def test_refused(): + mask = np.ones(SHAPE, np.uint8) + with pytest.raises(ValueError, match="signed integer"): + b.chunk2sparse(mask, dtype=np.uint16, negatives="list") + with pytest.raises(ValueError, match="signed integer"): + b.chunk2sparse(mask, dtype=np.float32, negatives="list") + with pytest.raises(ValueError, match="howmany"): + b.chunk2sparse(mask, dtype=np.int32, negatives="list", howmany=100) + with pytest.raises(ValueError, match='"drop" or "list"'): + b.chunk2sparse(mask, dtype=np.int32, negatives="keep") + assert b.chunk2sparse(mask, dtype=np.int32).negatives == "drop" From 36fb446a89d8433751c953616986ce40dd2c25ae Mon Sep 17 00:00:00 2001 From: Jonathan Wright Date: Sat, 14 Nov 2026 09:05:25 +0100 Subject: [PATCH 5/5] f32 pixels in the SIMD collects f32 pixels (dark-subtracted CCD frames, usually through from_frames) took the scalar collect. The six 32-bit kernels take a second template argument F: the same kernel on the same bits with a float compare (AVX-512 cmp_ps GT_OQ, SSE2/AVX2 cmpgt_ps without the sign bias, vcgtq_f32, the VSX float compares, the scalar tails through bslz4_cmp); NaN is never above the cut, as in the scalar loop. collect_gt calls them with F; collect_nz is the u32 one on the bits (every non-zero float, NaN included as before, and -0.0, which adds nothing to the dot and is never above a cut). f64 stays scalar. Zen 4, ceo2 frames as f32 minus a dark of 10, sparsify at cut 150 through from_frames: 355 -> 1002 fps (u32: 991). test_collect_tiers runs every value for f32 with negatives, -0.0, NaN and +-inf and fractional cuts. x86 1571 passed; POWER9 1409 passed; aarch64 cross-build clean. Co-Authored-By: Claude Opus 5.5 --- src/_pipeline.py | 18 ++--- src/pipeline/pipeline.h | 2 +- src/pipeline/registry.c | 2 +- src/steps/collect.hpp | 149 ++++++++++++++++++++++++++++--------- test/test_collect_tiers.py | 41 ++++++++-- 5 files changed, 159 insertions(+), 53 deletions(-) diff --git a/src/_pipeline.py b/src/_pipeline.py index eacacc7..fd6581b 100644 --- a/src/_pipeline.py +++ b/src/_pipeline.py @@ -56,13 +56,13 @@ ), "collect": ( "auto", - ("avx512cs", "u8/u16/u32/i8/i16/i32 pixels; AVX-512 BW/VL and VBMI2"), - ("avx2cs", "u8/u16/u32/i8/i16/i32 pixels; AVX2"), - ("vsx", "u8/u16/u32/i8/i16/i32 pixels; POWER VSX"), - ("neon", "u8/u16/u32/i8/i16/i32 pixels; aarch64 NEON"), - ("sse2", "u8/u16/u32/i8/i16/i32 pixels; SSE2"), + ("avx512cs", "u8/u16/u32/i8/i16/i32/f32 pixels; AVX-512 BW/VL and VBMI2"), + ("avx2cs", "u8/u16/u32/i8/i16/i32/f32 pixels; AVX2"), + ("vsx", "u8/u16/u32/i8/i16/i32/f32 pixels; POWER VSX"), + ("neon", "u8/u16/u32/i8/i16/i32/f32 pixels; aarch64 NEON"), + ("sse2", "u8/u16/u32/i8/i16/i32/f32 pixels; SSE2"), ("scalar", "any pixel type"), - ("avx512f", "u8/u16/u32/i8/i16/i32 pixels; AVX-512 F/BW, for CPUs without VBMI2"), + ("avx512f", "u8/u16/u32/i8/i16/i32/f32 pixels; AVX-512 F/BW, for CPUs without VBMI2"), ), "route": ( "auto", @@ -98,8 +98,8 @@ ("f", 4): "f32", ("f", 8): "f64"} -# the pixel dtypes of the SIMD collects (signed ones compare as signed) -SIMD_COLLECT_DTYPES = ("u8", "u16", "u32", "i8", "i16", "i32") +# the pixel dtypes of the SIMD collects (signed and float ones compare as such) +SIMD_COLLECT_DTYPES = ("u8", "u16", "u32", "i8", "i16", "i32", "f32") # the pixel dtypes of the low-planes untransposes (i16: a block with a negative pixel takes kcb) LOWPLANES_DTYPES = ("u8", "u16", "i16") @@ -242,7 +242,7 @@ def resolve(pipeline, dtype, codec, mask, dot_auto=None, matrix=False): % (uname, np.dtype(dtype))) cname = NAMES["collect"][vals[COLLECT]] if cname != "scalar" and suffix not in SIMD_COLLECT_DTYPES: - raise ValueError("pipeline['collect']=%r is for u8/u16/u32/i8/i16/i32 pixels, not %s (use scalar)" + raise ValueError("pipeline['collect']=%r is for u8/u16/u32/i8/i16/i32/f32 pixels, not %s (use scalar)" % (cname, np.dtype(dtype))) if matrix: if vals[DOT] == 0: diff --git a/src/pipeline/pipeline.h b/src/pipeline/pipeline.h index c1cb70b..97d07cf 100644 --- a/src/pipeline/pipeline.h +++ b/src/pipeline/pipeline.h @@ -75,7 +75,7 @@ enum { #define BSLZ4_IS_LOWPLANES(u) (((u) >= BSLZ4_UNTRANSPOSE_LOWPLANES_VBMI && (u) <= BSLZ4_UNTRANSPOSE_LOWPLANES_C) || \ (u) == BSLZ4_UNTRANSPOSE_LOWPLANES_AVX512BW) -/* collect: the SIMD values take u8/u16/u32 and i8/i16/i32 (a signed compare) */ +/* collect: the SIMD values take u8/u16/u32, i8/i16/i32 (a signed compare) and f32 (a float compare) */ /* negatives="list" (signed integer pixels): OR'ed into the pipeline's collect value by Python and * into bslz4_work.collect by the driver -- the pixel list also holds every negative pixel (for diff --git a/src/pipeline/registry.c b/src/pipeline/registry.c index 2f04b87..31e70fd 100644 --- a/src/pipeline/registry.c +++ b/src/pipeline/registry.c @@ -26,7 +26,7 @@ typedef char bslz4_value_slots_cover_every_step #define BSLZ4_ALL_DTYPES 0x3FFu #define BSLZ4_INT_DTYPES 0x0FFu /* u8 .. i64 */ #define BSLZ4_INT32_DTYPES 0x077u /* u8 u16 u32 i8 i16 i32: sums of 2**31 values fit int64 */ -#define BSLZ4_SIMD_COLLECT_DTYPES 0x077u /* the SIMD collects: u8 u16 u32 i8 i16 i32 (signed: a signed compare) */ +#define BSLZ4_SIMD_COLLECT_DTYPES 0x177u /* the SIMD collects: u8 u16 u32 i8 i16 i32 f32 (signed, float: their own compare) */ #define BSLZ4_LOWPLANES_DTYPES ((1u << 0) | (1u << 1) | (1u << 5)) /* the low-planes values: u8, and u16 / i16 blocks in 0..255 (a negative i16 sets the high planes) */ diff --git a/src/steps/collect.hpp b/src/steps/collect.hpp index 87e777f..2e8d6e4 100644 --- a/src/steps/collect.hpp +++ b/src/steps/collect.hpp @@ -5,8 +5,8 @@ * sparse-route compaction (mask>0 & value!=0, called here with cut==0 -- * identical to !=0 for unsigned T). * - * Six SIMD collect values, for u8/u16/u32 and i8/i16/i32 pixels (signed ones: the - * same kernels on the same bits with a signed compare, template argument S): avx512cs (AVX-512 F+BW+VL+VBMI2, + * Six SIMD collect values, for u8/u16/u32, i8/i16/i32 and f32 pixels (signed and float ones: + * the same kernels on the same bits with a signed (S) or float (F) compare): avx512cs (AVX-512 F+BW+VL+VBMI2, * compress-store), avx512f (AVX-512 F+BW+VL without VBMI2: dword compress), * avx2cs (AVX2, table-driven compaction), sse2 (x86-64 ABI * baseline, always present), VSX (POWER8+, probe in @@ -28,8 +28,8 @@ * BSLZ4_COLLECT_VSX POWER8+ (needs -maltivec -mvsx); u8/u16/u32. * BSLZ4_COLLECT_NEON AArch64 ASIMD (ARMv8-A baseline); u8/u16/u32. * BSLZ4_COLLECT_SSE2 x86-64 baseline, always present; u8/u16/u32. - * BSLZ4_COLLECT_SCALAR the plain scalar loop, any dtype; u64, i64 and the - * floats take it whatever the value. + * BSLZ4_COLLECT_SCALAR the plain scalar loop, any dtype; u64, i64 and f64 + * take it whatever the value. * * The value actually used is pipeline[BSLZ4_STEP_COLLECT], validated by * bslz4_resolve (src/pipeline/registry.c). Availability is decided by @@ -93,6 +93,23 @@ inline bool bslz4_gt(U a, U b) { return S ? (I) a > (I) b : a > b; } +/* v > cut on the bits of float pixels (F): the 32-bit kernels serve f32 too */ +template struct bslz4_fgt { + static inline bool gt(U, U) { return false; } +}; +template<> struct bslz4_fgt { + static inline bool gt(uint32_t a, uint32_t b) { + float x, y; + memcpy(&x, &a, 4); + memcpy(&y, &b, 4); + return x > y; + } +}; +template +inline bool bslz4_cmp(U a, U b) { + return F ? bslz4_fgt::gt(a, b) : bslz4_gt(a, b); +} + /* v < 0, false for unsigned and float-free of -Wtype-limits */ template::value> struct bslz4_neg { static inline bool of(T) { return false; } @@ -109,29 +126,31 @@ inline bool bslz4_keep(T v, T cut, int collect) { #if BSLZ4_HAVE_VSX_COLLECT /* the VSX compares, unsigned or (S) signed */ -template +template inline __vector __bool char bslz4_vsx_gt(__vector unsigned char v, __vector unsigned char c) { return S ? vec_cmpgt((__vector signed char) v, (__vector signed char) c) : vec_cmpgt(v, c); } -template +template inline int bslz4_vsx_any_gt(__vector unsigned char v, __vector unsigned char c) { return S ? vec_any_gt((__vector signed char) v, (__vector signed char) c) : vec_any_gt(v, c); } -template +template inline __vector __bool short bslz4_vsx_gt(__vector unsigned short v, __vector unsigned short c) { return S ? vec_cmpgt((__vector signed short) v, (__vector signed short) c) : vec_cmpgt(v, c); } -template +template inline int bslz4_vsx_any_gt(__vector unsigned short v, __vector unsigned short c) { return S ? vec_any_gt((__vector signed short) v, (__vector signed short) c) : vec_any_gt(v, c); } -template +template inline __vector __bool int bslz4_vsx_gt(__vector unsigned int v, __vector unsigned int c) { - return S ? vec_cmpgt((__vector signed int) v, (__vector signed int) c) : vec_cmpgt(v, c); + return F ? vec_cmpgt((__vector float) v, (__vector float) c) + : S ? vec_cmpgt((__vector signed int) v, (__vector signed int) c) : vec_cmpgt(v, c); } -template +template inline int bslz4_vsx_any_gt(__vector unsigned int v, __vector unsigned int c) { - return S ? vec_any_gt((__vector signed int) v, (__vector signed int) c) : vec_any_gt(v, c); + return F ? vec_any_gt((__vector float) v, (__vector float) c) + : S ? vec_any_gt((__vector signed int) v, (__vector signed int) c) : vec_any_gt(v, c); } #endif @@ -267,7 +286,7 @@ inline int bslz4_collect_avx512cs_u16(const uint16_t *BSLZ4_RESTRICT block, /* u32 version of collect avx512cs: vpcompressd is plain AVX-512F. Same store * bound as the u16 one, with 16 entries per group. */ -template +template __attribute__((target("avx512f,avx512bw,avx512vl,popcnt"))) inline int bslz4_collect_avx512cs_u32(const uint32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, @@ -280,7 +299,8 @@ inline int bslz4_collect_avx512cs_u32(const uint32_t *BSLZ4_RESTRICT block, const __m512i iota = _mm512_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15); for (; j + 16 <= block_elems; j += 16) { const __m512i v = _mm512_loadu_si512((const void *) &block[j]); - __mmask16 m = (S ? _mm512_cmpgt_epi32_mask(v, vcut) : _mm512_cmpgt_epu32_mask(v, vcut)); + __mmask16 m = (F ? _mm512_cmp_ps_mask(_mm512_castsi512_ps(v), _mm512_castsi512_ps(vcut), _CMP_GT_OQ) + : S ? _mm512_cmpgt_epi32_mask(v, vcut) : _mm512_cmpgt_epu32_mask(v, vcut)); if (mask) m &= _mm_cmpneq_epu8_mask(_mm_loadu_si128((const __m128i *) &mask[j + i0]), _mm_setzero_si128()); @@ -291,7 +311,7 @@ inline int bslz4_collect_avx512cs_u32(const uint32_t *BSLZ4_RESTRICT block, npx += __builtin_popcount((unsigned) m); } for (; j < block_elems; j++) { - if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_cmp(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -374,12 +394,12 @@ inline int bslz4_avx512f_emit16(__m512i v32, __mmask16 mq, size_t base, int npx, return npx + __builtin_popcount((unsigned) mq); } -template +template inline int bslz4_collect_tail(const T *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t j, size_t n, T cut, int npx, T *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { for (; j < n; j++) { - if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_cmp(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -457,7 +477,7 @@ inline int bslz4_collect_avx512f_u16(const uint16_t *BSLZ4_RESTRICT block, const return bslz4_collect_tail(block, mask, i0, j, n, cut, npx, out_vals, out_adr); } -template +template __attribute__((target("avx512f,avx512bw,popcnt"))) inline int bslz4_collect_avx512f_u32(const uint32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, uint32_t cut, @@ -467,14 +487,15 @@ inline int bslz4_collect_avx512f_u32(const uint32_t *BSLZ4_RESTRICT block, const const __m512i vcut = _mm512_set1_epi32((int) cut); for (; j + 16 <= n; j += 16) { const __m512i v = _mm512_loadu_si512((const void *) &block[j]); - __mmask16 m = (S ? _mm512_cmpgt_epi32_mask(v, vcut) : _mm512_cmpgt_epu32_mask(v, vcut)); + __mmask16 m = (F ? _mm512_cmp_ps_mask(_mm512_castsi512_ps(v), _mm512_castsi512_ps(vcut), _CMP_GT_OQ) + : S ? _mm512_cmpgt_epi32_mask(v, vcut) : _mm512_cmpgt_epu32_mask(v, vcut)); if (mask) m &= (__mmask16) ~(uint32_t) _mm_movemask_epi8( _mm_cmpeq_epi8(_mm_loadu_si128((const __m128i *) (const void *) &mask[j + i0]), _mm_setzero_si128())); if (!m) continue; npx = bslz4_avx512f_emit16(v, m, j + i0, npx, out_vals, out_adr); } - return bslz4_collect_tail(block, mask, i0, j, n, cut, npx, out_vals, out_adr); + return bslz4_collect_tail(block, mask, i0, j, n, cut, npx, out_vals, out_adr); } #endif /* BSLZ4_HAVE_AVX512_COLLECT */ @@ -593,7 +614,7 @@ inline int bslz4_collect_avx2cs_u16(const uint16_t *BSLZ4_RESTRICT block, return npx; } -template +template __attribute__((target("avx2"))) inline int bslz4_collect_avx2cs_u32(const uint32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, @@ -602,7 +623,7 @@ inline int bslz4_collect_avx2cs_u32(const uint32_t *BSLZ4_RESTRICT block, uint32_t *BSLZ4_RESTRICT out_adr) { int npx = 0; size_t j = 0; - const __m256i sign_bias = _mm256_set1_epi32((int) (S ? 0 : 0x80000000u)); + const __m256i sign_bias = _mm256_set1_epi32((int) ((S || F) ? 0 : 0x80000000u)); const __m256i vcut_b = _mm256_xor_si256(_mm256_set1_epi32((int) cut), sign_bias); for (; j + 8 <= block_elems; j += 8) { uint32_t mnz8 = mask @@ -611,7 +632,8 @@ inline int bslz4_collect_avx2cs_u32(const uint32_t *BSLZ4_RESTRICT block, __m256i v = _mm256_loadu_si256((const __m256i *) &block[j]); __m256i vb = _mm256_xor_si256(v, sign_bias); - __m256i gt = _mm256_cmpgt_epi32(vb, vcut_b); + __m256i gt = F ? _mm256_castps_si256(_mm256_cmp_ps(_mm256_castsi256_ps(vb), _mm256_castsi256_ps(vcut_b), _CMP_GT_OQ)) + : _mm256_cmpgt_epi32(vb, vcut_b); uint32_t mgt8 = (uint32_t) _mm256_movemask_ps(_mm256_castsi256_ps(gt)); uint32_t bits = mnz8 & mgt8; @@ -624,7 +646,7 @@ inline int bslz4_collect_avx2cs_u32(const uint32_t *BSLZ4_RESTRICT block, } } for (; j < block_elems; j++) { - if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_cmp(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -744,7 +766,7 @@ inline int bslz4_collect_sse2_u16(const uint16_t *BSLZ4_RESTRICT block, return npx; } -template +template __attribute__((target("sse2"))) inline int bslz4_collect_sse2_u32(const uint32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, @@ -753,7 +775,7 @@ inline int bslz4_collect_sse2_u32(const uint32_t *BSLZ4_RESTRICT block, uint32_t *BSLZ4_RESTRICT out_adr) { int npx = 0; size_t j = 0; - const __m128i sign_bias = _mm_set1_epi32((int) (S ? 0 : 0x80000000u)); + const __m128i sign_bias = _mm_set1_epi32((int) ((S || F) ? 0 : 0x80000000u)); const __m128i vcut_b = _mm_xor_si128(_mm_set1_epi32((int) cut), sign_bias); for (; j + 4 <= block_elems; j += 4) { uint32_t mnz4 = 0xFu; @@ -765,7 +787,8 @@ inline int bslz4_collect_sse2_u32(const uint32_t *BSLZ4_RESTRICT block, __m128i v = _mm_loadu_si128((const __m128i *) &block[j]); __m128i vb = _mm_xor_si128(v, sign_bias); - __m128i gt = _mm_cmpgt_epi32(vb, vcut_b); + __m128i gt = F ? _mm_castps_si128(_mm_cmpgt_ps(_mm_castsi128_ps(vb), _mm_castsi128_ps(vcut_b))) + : _mm_cmpgt_epi32(vb, vcut_b); uint32_t mgt4 = (uint32_t) _mm_movemask_ps(_mm_castsi128_ps(gt)); uint32_t bits = mnz4 & mgt4; @@ -778,7 +801,7 @@ inline int bslz4_collect_sse2_u32(const uint32_t *BSLZ4_RESTRICT block, } } for (; j < block_elems; j++) { - if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_cmp(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -921,7 +944,7 @@ inline int bslz4_collect_vsx_u16(const uint16_t *BSLZ4_RESTRICT block, return npx; } -template +template inline int bslz4_collect_vsx_u32(const uint32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, uint32_t cut, @@ -937,7 +960,7 @@ inline int bslz4_collect_vsx_u32(const uint32_t *BSLZ4_RESTRICT block, const int start = npx; if (dense) { for (; j < end; j++) { - if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_cmp(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -947,8 +970,8 @@ inline int bslz4_collect_vsx_u32(const uint32_t *BSLZ4_RESTRICT block, for (; j < end; j += 4) { __vector unsigned int v; memcpy(&v, &block[j], sizeof(v)); - if (!bslz4_vsx_any_gt(v, vcut)) continue; - __vector __bool int gt = bslz4_vsx_gt(v, vcut); + if (!bslz4_vsx_any_gt(v, vcut)) continue; + __vector __bool int gt = bslz4_vsx_gt(v, vcut); for (int lane = 0; lane < 4; lane++) { if (bslz4_mask_ok(mask, j + i0 + lane) && vec_extract((__vector unsigned int) gt, lane)) { out_vals[npx] = vec_extract(v, lane); @@ -961,7 +984,7 @@ inline int bslz4_collect_vsx_u32(const uint32_t *BSLZ4_RESTRICT block, dense = bslz4_vsx_next_dense(dense, npx - start); } for (; j < n; j++) { - if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_cmp(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -1075,7 +1098,7 @@ inline int bslz4_collect_neon_u16(const uint16_t *BSLZ4_RESTRICT block, return npx; } -template +template inline int bslz4_collect_neon_u32(const uint32_t *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, size_t i0, size_t n, uint32_t cut, @@ -1086,7 +1109,8 @@ inline int bslz4_collect_neon_u32(const uint32_t *BSLZ4_RESTRICT block, uint32x4_t vcut = vdupq_n_u32(cut); for (; j + 4 <= n; j += 4) { uint32x4_t v = vld1q_u32(&block[j]); - uint32x4_t gt = (S ? vcgtq_s32(vreinterpretq_s32_u32(v), vreinterpretq_s32_u32(vcut)) : vcgtq_u32(v, vcut)); + uint32x4_t gt = (F ? vcgtq_f32(vreinterpretq_f32_u32(v), vreinterpretq_f32_u32(vcut)) + : S ? vcgtq_s32(vreinterpretq_s32_u32(v), vreinterpretq_s32_u32(vcut)) : vcgtq_u32(v, vcut)); if (vmaxvq_u32(gt) == 0) continue; uint32_t vbuf[4], gbuf[4]; vst1q_u32(vbuf, v); @@ -1100,7 +1124,7 @@ inline int bslz4_collect_neon_u32(const uint32_t *BSLZ4_RESTRICT block, } } for (; j < n; j++) { - if (bslz4_mask_ok(mask, j + i0) & bslz4_gt(block[j], cut)) { + if (bslz4_mask_ok(mask, j + i0) & bslz4_cmp(block[j], cut)) { out_vals[npx] = block[j]; out_adr[npx] = (uint32_t) (j + i0); npx++; @@ -1577,6 +1601,59 @@ inline int bslz4_collect_nz(const int32_t *BSLZ4_RESTRICT block, const (uint32_t *) (void *) out_vals, out_adr); } +/* f32 pixels: the 32-bit SIMD collects with a float compare */ + +template<> +inline int bslz4_collect_gt(const float *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, + size_t i0, size_t n, float cut, int collect, + float *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { + /* the 32-bit kernels on the same bits, comparing as floats (NaN is never above the cut) */ + const uint32_t *ub = (const uint32_t *) (const void *) block; + uint32_t *uo = (uint32_t *) (void *) out_vals; + uint32_t uc; + memcpy(&uc, &cut, 4); + (void) ub; (void) uo; (void) uc; +#if BSLZ4_HAVE_AVX512_COLLECT +#if BSLZ4_HAVE_VBMI_GFNI + if (collect == BSLZ4_COLLECT_AVX512CS) return bslz4_collect_avx512cs_u32(ub, mask, i0, n, uc, uo, out_adr); +#endif +#endif +#if BSLZ4_HAVE_AVX512_COLLECT + if (collect == BSLZ4_COLLECT_AVX512F) return bslz4_collect_avx512f_u32(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_AVX2_COLLECT + if (collect == BSLZ4_COLLECT_AVX2CS) return bslz4_collect_avx2cs_u32(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_SSE2_COLLECT + if (collect == BSLZ4_COLLECT_SSE2) return bslz4_collect_sse2_u32(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_VSX_COLLECT + if (collect == BSLZ4_COLLECT_VSX) return bslz4_collect_vsx_u32(ub, mask, i0, n, uc, uo, out_adr); +#endif +#if BSLZ4_HAVE_NEON_COLLECT + if (collect == BSLZ4_COLLECT_NEON) return bslz4_collect_neon_u32(ub, mask, i0, n, uc, uo, out_adr); +#endif + int npx = 0; + for (size_t j = 0; j < n; j++) { + if (BSLZ4_UNLIKELY(bslz4_mask_ok(mask, j + i0) & (block[j] > cut))) { + out_vals[npx] = block[j]; + out_adr[npx] = (uint32_t) (j + i0); + npx++; + } + } + return npx; +} + +template<> +inline int bslz4_collect_nz(const float *BSLZ4_RESTRICT block, const uint8_t *BSLZ4_RESTRICT mask, + size_t i0, size_t n, int collect, + float *BSLZ4_RESTRICT out_vals, uint32_t *BSLZ4_RESTRICT out_adr) { + /* the bits != 0: every non-zero float, NaN included as in the scalar loop, and -0.0, which + * adds nothing to the sparse route's dot and is never above the cut */ + return bslz4_collect_nz((const uint32_t *) (const void *) block, mask, i0, n, collect, + (uint32_t *) (void *) out_vals, out_adr); +} + #endif /* any tier available */ } /* namespace bslz4 */ diff --git a/test/test_collect_tiers.py b/test/test_collect_tiers.py index 9387b5b..9f7fd68 100644 --- a/test/test_collect_tiers.py +++ b/test/test_collect_tiers.py @@ -33,8 +33,10 @@ if i and _pipeline.available("collect", i) == 1] -def make(kind, dt, rng): +def make(kind, dt, rng, specials=True): n = (3,) + SHAPE + if np.dtype(dt).kind == "f": + return make_float(kind, dt, rng, specials) top = np.iinfo(dt).max if kind == "sparse": a = np.zeros(n, np.int64) @@ -53,7 +55,30 @@ def make(kind, dt, rng): return a.astype(dt) -DTYPES = [np.uint8, np.uint16, np.uint32, np.int8, np.int16, np.int32] +def make_float(kind, dt, rng, specials): + """Fractional values with negatives (dark-subtracted frames); specials: also -0.0, NaN and + +-inf, which numpy and the kernels both never count as above a finite cut (but inf).""" + n = (3,) + SHAPE + if kind == "sparse": + a = np.zeros(n) + a[rng.random(n) < 1e-3] = 900.5 + a[:, 10:14, 50:54] = 900.25 + elif kind == "medium": + a = rng.poisson(0.1, n) + rng.normal(0, 0.3, n) + elif kind == "dense": + a = rng.normal(30, 10, n) + else: # "max" + a = np.full(n, 3.0e38) + a = a.astype(dt) + if specials: + a[rng.random(n) < 0.005] = -0.0 + a[rng.random(n) < 0.002] = np.nan + a[rng.random(n) < 0.001] = np.inf + a[rng.random(n) < 0.001] = -np.inf + return a + + +DTYPES = [np.uint8, np.uint16, np.uint32, np.int8, np.int16, np.int32, np.float32] def chunks_of(frames, path): @@ -72,8 +97,12 @@ def test_every_tier_matches_numpy(tmp_path, dt, kind, masked): chunks, codec = chunks_of(frames, tmp_path / "c.h5") mask = ((rng.random(SHAPE) > 0.15) if masked else np.ones(SHAPE, bool)).astype(np.uint8) # the threshold is a C int: the largest cut tried is min(max - 1, 2**31 - 1) - top = min(int(np.iinfo(dt).max) - 1, 2 ** 31 - 1) - for cut in sorted(set(min(c, top) for c in (0, 3, 600, top))): + if np.dtype(dt).kind == "f": + cuts = (0, 0.25, 3, 600.5, 1e30) + else: + top = min(int(np.iinfo(dt).max) - 1, 2 ** 31 - 1) + cuts = sorted(set(min(c, top) for c in (0, 3, 600, top))) + for cut in cuts: for tier in TIERS: p = {"collect": tier, "untranspose": "kcb", "mask": "pixel" if masked else "none"} @@ -85,14 +114,14 @@ def test_every_tier_matches_numpy(tmp_path, dt, kind, masked): n = int(npx[f]) assert n == want.size, (tier, kind, cut, f) assert (adr[f, :n] == want).all(), (tier, kind, cut, f) - assert (vals[f, :n] == flat[want]).all(), (tier, kind, cut, f) + assert np.array_equal(vals[f, :n], flat[want]), (tier, kind, cut, f) @pytest.mark.parametrize("dt", DTYPES) @pytest.mark.parametrize("route", ["dense", "sparse"]) def test_powder_same_for_every_tier(tmp_path, dt, route): rng = np.random.default_rng(22) - frames = make("medium", dt, rng) + frames = make("medium", dt, rng, specials=False) chunks, codec = chunks_of(frames, tmp_path / "p.h5") npix = SHAPE[0] * SHAPE[1] mask = (rng.random(SHAPE) > 0.1).astype(np.uint8)