diff --git a/docs/index.md b/docs/index.md index 3d3fc41..37a6e24 100644 --- a/docs/index.md +++ b/docs/index.md @@ -35,6 +35,34 @@ 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) 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 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..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,10 +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. + 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 @@ -516,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 @@ -584,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 @@ -728,8 +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 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 @@ -751,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): @@ -879,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) @@ -997,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.py b/src/_pipeline.py index 1bfe1e6..fd6581b 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/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 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,6 +98,12 @@ ("f", 4): "f32", ("f", 8): "f64"} +# 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") + + 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/f32 pixels, not %s (use scalar)" % (cname, np.dtype(dtype))) if matrix: if vals[DOT] == 0: 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 71a58d7..97d07cf 100644 --- a/src/pipeline/pipeline.h +++ b/src/pipeline/pipeline.h @@ -75,7 +75,13 @@ 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, 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 + * 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 614196c..31e70fd 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 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) */ /* 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 @@ -164,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; @@ -177,9 +180,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]]; @@ -188,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 11d109a..2e8d6e4 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, 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 @@ -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 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 @@ -50,6 +51,7 @@ #include "../pipeline/pipeline.h" #include +#include #include "c2py_amd64.h" #include "c2py_arm64.h" @@ -82,6 +84,76 @@ 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; +} + +/* 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; } +}; +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 +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 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 +inline int bslz4_vsx_any_gt(__vector unsigned int v, __vector unsigned int 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 + /* 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 +246,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 +260,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 +275,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 +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 __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 +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 = _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()); @@ -236,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) & (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++; @@ -250,6 +325,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 +339,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 +355,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 +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) & (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++; @@ -332,6 +408,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 +418,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 +442,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 +455,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 +474,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 +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 = _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 */ @@ -476,6 +556,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 +566,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 +605,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 +614,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 +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) 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 @@ -550,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; @@ -563,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) & (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++; @@ -575,6 +658,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 +668,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 +702,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 +726,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 +735,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 +757,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 +766,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 +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) 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; @@ -701,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; @@ -714,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) & (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++; @@ -724,6 +811,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 +820,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 +838,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 +895,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 +911,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 +921,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 +935,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 +944,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 +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) & (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++; @@ -880,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 (!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 +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) & (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++; @@ -904,6 +994,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 +1010,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 +1020,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 +1034,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 +1064,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 +1075,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 +1089,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 +1098,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 +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 = 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); @@ -1030,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) & (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++; @@ -1040,6 +1134,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 +1145,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 +1159,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++; @@ -1087,10 +1182,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++; @@ -1342,6 +1437,223 @@ 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; + /* 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 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 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 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 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 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 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) | (neg & (block[j] < 0))))) { + 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 & ~BSLZ4_COLLECT_NEGATIVES, + (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; + /* 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 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 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 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 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 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 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) | (neg & (block[j] < 0))))) { + 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 & ~BSLZ4_COLLECT_NEGATIVES, + (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; + /* 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 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 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 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 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 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 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) | (neg & (block[j] < 0))))) { + 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 & ~BSLZ4_COLLECT_NEGATIVES, + (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/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_collect_tiers.py b/test/test_collect_tiers.py index 59a5b38..9f7fd68 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. @@ -32,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) @@ -45,7 +48,37 @@ 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) + + +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): @@ -55,7 +88,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): @@ -64,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"} @@ -77,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", [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) - 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) 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"}) 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"