Cut peak memory across the encode and payload paths

CHANGELOG.md

@@ -36,6 +36,50 @@ in the README).
   same validated style properties, so an explicit `style=` still wins and specs
   that don't use them are byte-identical.
 
+### Changed
+- Peak memory cut on four paths, with byte-identical output everywhere (89
+  payload/export/view fingerprints pinned before and after):
+  - The indexed-palette PNG encoder stages one scanline buffer and hands it to
+    zlib directly instead of building a per-row `bytes` list, joining it, and
+    narrowing an `intp`-per-pixel `np.unique` inverse. A 1800×840 export peaks
+    at 6 MB instead of 50 MB (33 MB instead of 274 MB at 4K) and encodes ~1.8x
+    faster; the truecolor branch drops ~17%.
+  - Full-column color quantization (density mean-color planes, `direct_rgba`
+    channels, u8 live-wire channels) runs chunk-bounded and in place. Its chunk
+    was 4M rows, so every real column still paid the one-shot peak: a 2.1M-row
+    continuous color channel resolved in 44 MB and now resolves in 7 MB, taking
+    a colored 2.1M-point first paint from 181 MB to 147 MB of RSS.
+  - Direct-tier scatter/line/area payloads skip the all-visible row mask. Zone
+    maps already count NaN *and* ±inf as null, so on linear axes with no nulls
+    the mask is provably all-true; it was three O(N) passes and two N-byte
+    temporaries per build. Emit is 16–35% faster, and area no longer allocates
+    an identity index vector (nor gathers animation keys through it).
+  - SVG documents assemble in one flat join with block-buffered markers, and the
+    native rasterizer borrows the display list instead of freezing a `bytes`
+    copy. A 100k-point SVG export peaks at 27 MB instead of 39 MB.
+  - The JPEG encoder streams instead of exploding: the entropy packer works in
+    bounded bit passes (it cost 17 bytes per output *bit*, so a 2.8 MB stream
+    peaked over 1.5 GB), the YCbCr planes are released as they are consumed, the
+    quantize chain rounds in place via `trunc(x + copysign(0.5, x))`, and the
+    per-component token fields are freed as they are gathered. A 1800x840 export
+    peaks at 48 MB instead of 108 (photographic: **400 MB instead of 1516 MB**)
+    and is 5-19% faster.
+  - Standalone HTML export joins the document once from parts, with every large
+    string — the client bundle, the spec, each base64 chunk — as its own part, so
+    the join copies it exactly once. Previously the chunks were folded through a
+    `"\n".join(...)` and then into an f-string, duplicating 4/3 of the payload. A
+    1M-point export peaks at 33 MB instead of 41; a small export (where the
+    ~330 KB client bundle is the document) is ~2% faster.
+  - The native rasterizer's serial point and segment passes no longer
+    materialize an `0..n` index vector to read back in order (4 bytes per mark),
+    and the per-point density quantizer keeps one temporary instead of three.
+- `memory_report()` reports `capacity_bytes` per column and
+  `canonical_capacity_bytes` per store, and builds `resident_array_bytes` from
+  the capacity total. A streamed column's `values` is a prefix view of its
+  capacity-doubling buffer, so up to half of what it holds was invisible to the
+  report (§27: if a number isn't in the report, it isn't real). Figures that
+  never appended report exactly what they did before.
+
 ### Fixed
 - The colorbar stringified its colormap, so a custom ramp reached it as an
   unparseable name and silently painted viridis while the marks beside it

python/xy/_figure.py

@@ -1871,8 +1871,10 @@ def memory_report(self) -> dict[str, Any]:
         report["transport_bytes_per_point"] = len(blob) / n_total
         report["pyramid_bytes"] = interaction.pyramid_report_bytes(self)
         report["bin_color_bytes"] = interaction.bin_color_cache_bytes(self)
+        # Capacity, not live length: a streamed column's growth-buffer slack is
+        # resident RAM (§27), and equals `canonical_bytes` when nothing appended.
         report["resident_array_bytes"] = (
-            report["canonical_bytes"]
+            report["canonical_capacity_bytes"]
             + report["channel_bytes"]
             + report["pyramid_bytes"]
             + report["bin_color_bytes"]

python/xy/_jpeg.py

@@ -20,7 +20,9 @@
 
 from __future__ import annotations
 
+import itertools
 import struct
+from typing import Optional
 
 import numpy as np
 
@@ -168,6 +170,9 @@ def _dct_matrix() -> np.ndarray:
 # searchsorted keeps this exact where a float log2 could round at the edges.
 _POW2 = np.int64(1) << np.arange(12, dtype=np.int64)
 
+#: Placeholder that releases a token field's storage while its list slot lives on.
+_EMPTY_I64 = np.empty(0, dtype=np.int64)
+
 
 def _bit_size(magnitude: np.ndarray) -> np.ndarray:
     return np.searchsorted(_POW2, magnitude, side="right").astype(np.int64)
@@ -255,22 +260,61 @@ def _component_tokens(
     return merged[0], merged[1], merged[2], merged[3], merged[4], merged[5]
 
 
+#: Bits packed per `_pack_entropy` pass. The exploded form costs 17 bytes per
+#: output *bit* (two int64 index vectors plus the bit itself), so packing a whole
+#: image at once peaked at ~130 MB for a 1 MB entropy stream. Bit packing is not
+#: chunk-independent — group boundaries land mid-byte — so each pass carries its
+#: leftover bits into the next one, which keeps the byte stream identical.
+_ENTROPY_BIT_CHUNK = 1 << 18
+
+
 def _pack_entropy(chunk: np.ndarray, nbits: np.ndarray) -> bytes:
     """MSB-first bit packing of (value, bit-count) chunks, with 1-padding to a
-    byte boundary and 0x00 stuffing after every 0xFF (both per T.81)."""
+    byte boundary and 0x00 stuffing after every 0xFF (both per T.81).
+
+    Packed in bounded passes: token groups of at most `_ENTROPY_BIT_CHUNK` bits,
+    with sub-byte remainders carried across group boundaries so the emitted bytes
+    are exactly those of a single whole-image pass. Byte stuffing is per-byte
+    local, so it applies per group.
+    """
     total = int(nbits.sum())
-    start = np.cumsum(nbits) - nbits
-    idx = np.repeat(np.arange(chunk.size, dtype=np.int64), nbits)
-    offset = np.arange(total, dtype=np.int64) - np.repeat(start, nbits)
-    bits = ((chunk[idx] >> (nbits[idx] - 1 - offset)) & 1).astype(np.uint8)
-    pad = (-total) % 8
-    if pad:
-        bits = np.concatenate((bits, np.ones(pad, dtype=np.uint8)))
-    stream = np.packbits(bits)
-    ff = np.flatnonzero(stream == 0xFF)
-    if ff.size:
-        stream = np.insert(stream, ff + 1, np.uint8(0))
-    return stream.tobytes()
+    cum = np.cumsum(nbits)
+    # Token index at each bit-budget boundary: a group is [t0, t1) tokens, and a
+    # single token never spans groups (its bits stay contiguous, as before).
+    edges = np.searchsorted(
+        cum, np.arange(_ENTROPY_BIT_CHUNK, total, _ENTROPY_BIT_CHUNK), side="left"
+    )
+    bounds = [0, *(int(e) + 1 for e in edges), chunk.size]
+    parts: list[bytes] = []
+    carry = np.empty(0, dtype=np.uint8)
+    for t0, t1 in itertools.pairwise(bounds):
+        if t1 <= t0:
+            continue
+        values = chunk[t0:t1]
+        widths = nbits[t0:t1]
+        group_total = int(widths.sum())
+        start = np.cumsum(widths) - widths
+        idx = np.repeat(np.arange(values.size, dtype=np.int64), widths)
+        offset = np.arange(group_total, dtype=np.int64) - np.repeat(start, widths)
+        bits = ((values[idx] >> (widths[idx] - 1 - offset)) & 1).astype(np.uint8)
+        if carry.size:
+            bits = np.concatenate((carry, bits))
+        whole = (bits.size // 8) * 8
+        if t1 == chunk.size:
+            # Final group: 1-pad the tail to a byte boundary (T.81).
+            pad = (-bits.size) % 8
+            if pad:
+                bits = np.concatenate((bits, np.ones(pad, dtype=np.uint8)))
+            whole = bits.size
+            carry = np.empty(0, dtype=np.uint8)
+        else:
+            carry = bits[whole:].copy()
+        stream = np.packbits(bits[:whole])
+        ff = np.flatnonzero(stream == 0xFF)
+        if ff.size:
+            stream = np.insert(stream, ff + 1, np.uint8(0))
+        parts.append(stream.tobytes())
+    return b"".join(parts)
 
 
 def _headers(h: int, w: int, qy_zz: np.ndarray, qc_zz: np.ndarray) -> bytes:
@@ -310,6 +354,26 @@ def _headers(h: int, w: int, qy_zz: np.ndarray, qc_zz: np.ndarray) -> bytes:
     return b"\xff\xd8" + app0 + dqt + sof0 + dht + sos
 
 
+def _gathered(fields: list[list[np.ndarray]], index: int, order: np.ndarray) -> np.ndarray:
+    """One token field, concatenated across components into MCU order.
+
+    The per-component pieces are released as they are joined, so a field costs
+    one joined copy plus the gathered result rather than both plus the parts —
+    each is an int64 per token, and a photographic image has millions. `parts`
+    is dropped explicitly rather than left to the frame: a local stays alive
+    until the function returns, so `return np.concatenate(parts)[order]` would
+    hold the pieces, the join, *and* the gather at once.
+    """
+    parts = [f[index] for f in fields]
+    for f in fields:
+        f[index] = _EMPTY_I64
+    joined = np.concatenate(parts)
+    del parts
+    gathered = joined[order]
+    del joined
+    return gathered
+
+
 def encode(rgba: np.ndarray, *, quality: int = 90) -> bytes:
     """Encode an `(h, w, 4)` RGBA (alpha ignored) or `(h, w, 3)` RGB uint8
     image as a baseline JFIF JPEG. Deterministic for identical input."""
@@ -338,6 +402,10 @@ def encode(rgba: np.ndarray, *, quality: int = 90) -> bytes:
     y = 0.299 * r + 0.587 * g + 0.114 * b - 128.0
     cb = -0.168736 * r - 0.331264 * g + 0.5 * b
     cr = 0.5 * r - 0.418688 * g - 0.081312 * b
+    # The interleaved f32 source is dead once the planes exist, and it is three
+    # floats per pixel — hand it back before the per-component pipeline runs.
+    # (`r`/`g`/`b` are views into it, so all four references have to go.)
+    del rgb, r, g, b
 
     qy = _scaled_quant(_QUANT_LUMA, quality)
     qc = _scaled_quant(_QUANT_CHROMA, quality)
@@ -346,37 +414,77 @@ def encode(rgba: np.ndarray, *, quality: int = 90) -> bytes:
 
     pad_h, pad_w = (-h) % 8, (-w) % 8
     h8, w8 = h + pad_h, w + pad_w
-    fields: list[tuple[np.ndarray, ...]] = []
+    fields: list[list[np.ndarray]] = []
     keys: list[np.ndarray] = []
-    for comp, (plane, q_zz) in enumerate(zip((y, cb, cr), (qy_zz, qc_zz, qc_zz), strict=True)):
+    # Planes are released as they are consumed (one float per pixel each), so the
+    # pipeline holds one component's working set instead of all three.
+    planes: list[Optional[np.ndarray]] = [y, cb, cr]
+    del y, cb, cr
+    for comp, q_zz in enumerate((qy_zz, qc_zz, qc_zz)):
+        plane = planes[comp]
+        planes[comp] = None
+        assert plane is not None
         # Edge replication avoids the ringing a zero/black pad would inject
         # into every border block.
         padded = np.pad(plane, ((0, pad_h), (0, pad_w)), mode="edge")
+        del plane
         blocks = padded.reshape(h8 // 8, 8, w8 // 8, 8).swapaxes(1, 2).reshape(-1, 8, 8)
+        del padded
         coef = _DCT @ blocks @ _DCT.T
-        scaled = coef.reshape(-1, 64)[:, _ZIGZAG] / q_zz.astype(np.float32)
+        del blocks
+        # The zigzag gather is a fresh array, so the divide, the round, and the
+        # sign restore all run in place on it: same f32 operations in the same
+        # order as the one-shot chain, five fewer full-size temporaries.
+        scaled = coef.reshape(-1, 64)[:, _ZIGZAG]
+        del coef
+        scaled /= q_zz.astype(np.float32)
         # Round half away from zero: any deterministic tie rule is valid JPEG;
-        # this one matches the common integer implementations.
-        quant = (np.sign(scaled) * np.floor(np.abs(scaled) + 0.5)).astype(np.int64)
+        # this one matches the common integer implementations. `trunc(x + ±0.5)`
+        # is that rule in three in-place passes and one temporary, where
+        # `sign(x) * floor(|x| + 0.5)` took four passes and four full-size
+        # temporaries. It is exact, not approximate: IEEE addition is
+        # sign-symmetric, so for x < 0 the sum is the exact negation of
+        # `|x| + 0.5`, and truncation toward zero then matches the floor of that
+        # magnitude. The ±0 cases land on ±0.0 either way and integer-cast to 0.
+        bias = np.copysign(np.float32(0.5), scaled)
+        scaled += bias
+        del bias
+        np.trunc(scaled, out=scaled)
+        quant = scaled.astype(np.int64)
+        del scaled
         # Exact-math AC magnitudes cap at 1020 (category 10); clamp is one-LSB
         # insurance against float rounding ever minting category 11, which the
         # baseline AC tables cannot code.
         quant[:, 1:] = np.clip(quant[:, 1:], -1023, 1023)
         dc_tbl, ac_tbl = (0, 1) if comp == 0 else (2, 3)
         toks = _component_tokens(quant, dc_tbl, ac_tbl)
+        del quant
         # 4:4:4 → one block per component per MCU, so the MCU-interleaved
         # order Y, Cb, Cr is a stable sort on (block, component, in-block seq).
         keys.append(((toks[0] * 3 + comp) << 8) | toks[1])
-        fields.append(toks)
+        # Only the four emitted fields outlive the key: the block and in-block
+        # sequence numbers are folded into it, and on a photographic image each
+        # dropped field is an int64 per token.
+        fields.append(list(toks[2:]))
+        del toks
 
     order = np.argsort(np.concatenate(keys), kind="stable")
-    tbl = np.concatenate([f[2] for f in fields])[order]
-    sym = np.concatenate([f[3] for f in fields])[order]
-    ampl = np.concatenate([f[4] for f in fields])[order]
-    abits = np.concatenate([f[5] for f in fields])[order]
+    keys.clear()
+    tbl = _gathered(fields, 0, order)
+    sym = _gathered(fields, 1, order)
+    ampl = _gathered(fields, 2, order)
+    abits = _gathered(fields, 3, order)
+    del fields
     code = _HUFF_CODES[tbl, sym]
     clen = _HUFF_LENS[tbl, sym]
-    # Huffman code then amplitude bits, as one ≤27-bit chunk per token.
-    entropy = _pack_entropy((code << abits) | ampl, clen + abits)
+    del tbl, sym
+    # Huffman code then amplitude bits, as one ≤27-bit chunk per token. Both
+    # gathers are fresh arrays, so the combine runs in place.
+    code <<= abits
+    code |= ampl
+    clen += abits
+    del ampl, abits
+    entropy = _pack_entropy(code, clen)
+    del code, clen
 
     return _headers(h, w, qy_zz, qc_zz) + entropy + b"\xff\xd9"

python/xy/_native.py

@@ -17,7 +17,7 @@
 import os
 import sys
 from pathlib import Path
-from typing import Optional
+from typing import Any, Optional
 
 import numpy as np
 import numpy.typing as npt
@@ -3093,8 +3093,11 @@ def _byte_span_arrays(spans):  # noqa: ANN001, ANN202 - private ctypes adapter
     return arenas, pointers, lengths
 
 
-def rasterize_spans(cmds: bytes, spans, w: int, h: int) -> npt.NDArray[np.uint8]:  # noqa: ANN001
-    """Paint a display list borrowing multiple call-scoped byte arenas."""
+def rasterize_spans(cmds: Any, spans, w: int, h: int) -> npt.NDArray[np.uint8]:  # noqa: ANN001
+    """Paint a display list borrowing multiple call-scoped byte arenas.
+
+    `cmds` is any read-only-safe buffer (`bytes`, `bytearray`, `memoryview`):
+    it is borrowed through `np.frombuffer`, never copied."""
     w = _positive_int(w, "raster width")
     h = _positive_int(h, "raster height")
     buf = np.frombuffer(cmds, dtype=np.uint8)
@@ -3115,8 +3118,9 @@ def rasterize_spans(cmds: bytes, spans, w: int, h: int) -> npt.NDArray[np.uint8]
     return out
 
 
-def rasterize_png_spans(cmds: bytes, spans, w: int, h: int) -> bytes:  # noqa: ANN001
-    """Paint and encode a display list borrowing multiple byte arenas."""
+def rasterize_png_spans(cmds: Any, spans, w: int, h: int) -> bytes:  # noqa: ANN001
+    """Paint and encode a display list borrowing multiple byte arenas
+    (`cmds` is borrowed, as in `rasterize_spans`)."""
     w = _positive_int(w, "raster width")
     h = _positive_int(h, "raster height")
     buf = np.frombuffer(cmds, dtype=np.uint8)

python/xy/_payload.py

@@ -403,13 +403,44 @@ def _finite_sel(t: Trace, xv: np.ndarray, yv: np.ndarray) -> np.ndarray | None:
             return None
         return np.flatnonzero(np.isfinite(xv) & np.isfinite(yv))
 
+    def _visible_mask_needed(
+        self,
+        t: Trace,
+        *,
+        prefiltered: bool,
+        base_column: Optional[Column] = None,
+    ) -> bool:
+        """Whether `_log_visible_mask` can drop any row for this trace.
+
+        The mask is three O(N) passes plus two N-byte temporaries, and on the
+        common shape (linear axes, no nulls) it is provably all-True: zone maps
+        count NaN *and* ±inf as null (§22/§19), so a null-free column has no
+        row for `isfinite` to reject, and `prefiltered` rows already went
+        through `_finite_sel`. Only a log axis (which additionally rejects
+        non-positive values) or a baseline column outside the x/y zone maps can
+        actually remove something.
+        """
+        if self._axis_scale(t.x_axis) == "log" or self._axis_scale(t.y_axis) == "log":
+            return True
+        if base_column is not None and base_column.zone.null_count:
+            return True
+        return not prefiltered and bool(t.x.zone.null_count or t.y.zone.null_count)
+
     def _log_visible_mask(
         self,
         t: Trace,
         xv: np.ndarray,
         yv: np.ndarray,
         base: Optional[np.ndarray] = None,
     ) -> np.ndarray:
+        """Rows this trace may actually ship.
+
+        Paired with `_visible_mask_needed`, which decides whether calling this
+        can drop anything at all: **a new rejection rule here needs a matching
+        condition there**, or the emitters will skip the mask on data it would
+        now reject. (`tests/test_figure.py` pins one case per rule; the
+        predicate going stale is otherwise silent.)
+        """
         mask = np.isfinite(xv) & np.isfinite(yv)
         if self._axis_scale(t.x_axis) == "log":
             mask &= xv > 0
@@ -483,7 +514,7 @@ def _emit_line(
             sel = self._finite_sel(t, xv, yv)
             if sel is not None:
                 xv, yv = xv[sel], yv[sel]
-        if len(xv):
+        if len(xv) and self._visible_mask_needed(t, prefiltered=sel is not None):
             finite = self._log_visible_mask(t, xv, yv)
             if not bool(np.all(finite)):
                 sel = np.flatnonzero(finite) if sel is None else sel[finite]
@@ -507,9 +538,11 @@ def _emit_area(
         tier, (xv, yv, bv) = self._m4_decimate(
             t, xr, px_width, t.x.values, t.y.values, t.base.values
         )
-        sel = np.flatnonzero(self._log_visible_mask(t, xv, yv, bv))
-        if len(sel) != len(xv):
-            xv, yv, bv = xv[sel], yv[sel], bv[sel]
+        sel = None
+        if self._visible_mask_needed(t, prefiltered=False, base_column=t.base):
+            sel = np.flatnonzero(self._log_visible_mask(t, xv, yv, bv))
+            if len(sel) != len(xv):
+                xv, yv, bv = xv[sel], yv[sel], bv[sel]
         entry = self._base_entry(t, pw, xv, yv, tier, self._default_styled(t))
         if tier == "decimated":
             entry["decimation_px"] = int(px_width)
@@ -535,7 +568,7 @@ def _emit_scatter(
         sel = self._finite_sel(t, xv, yv)
         if sel is not None:
             xv, yv = xv[sel], yv[sel]
-        if len(xv):
+        if len(xv) and self._visible_mask_needed(t, prefiltered=sel is not None):
             visible = self._log_visible_mask(t, xv, yv)
             if not bool(np.all(visible)):
                 sel = np.flatnonzero(visible) if sel is None else sel[visible]