Kernel gaps

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

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)

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:

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;

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

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;
 }