Performance: vectorizing generate_insar_mask() makes prepare_insar_hdf5 5x faster

[s-sasaki-earthsea-wizard]

Summary

generate_insar_mask() builds the InSAR mask layer with a pure-Python double loop over every output pixel — two scalar pybind SubSwaths.get_sample_sub_swath calls, NumPy scalar indexing and a list.append per pixel, ~5 µs/px. The per-pixel work is elementwise with no loop-carried dependency, so the same computation vectorizes with NumPy using only the existing SubSwaths.get_valid_samples_array bulk API. On a frame-scale grid (72.4 Mpx) the function drops from ~383 s to ~7 s on purely synthetic data (57x; self-contained reproducer below), and the containing prepare_insar_hdf5 stage runs 5x faster with bitwise-identical output on a real production run. The rewrite also simplifies the function body (guard clauses and a named membership helper instead of five levels of nesting). A candidate implementation is linked below; I can submit a PR with focused regression tests.

Mechanism

The inner loop (utils.py#L561-L616) executes per output pixel:

• two scalar pybind calls to SubSwaths.get_sample_sub_swath (reference, and offset-shifted secondary) via _compute_subswath_mask_id;
• NumPy scalar indexing into the two inputDataExceptionMask arrays plus Python-level bit packing;
• a list.append.

In a profiled production-scale run (details under Measurement), py-spy attributes 74.0% of all prepare_insar_hdf5 samples to the largest single mask call (the 6840×10581 RIFG interferogram grid) and 89.0% to generate_insar_mask overall (RUNW igram 11.3%, pixelOffsets masks 3.7%). The observed cost broadly tracks pixel count (RIFG/RUNW sample ratio 6.55 vs pixel ratio 6.94). Because the dominant term is interpreter-bound Python, the stage shows CPU/GPU parity (0.99x), and its dominant cost is non-GPU-addressable. Full profile (raw stacks, attribution tables, journal logs): profiling evidence.

Reproducing on synthetic data

No NISAR granules are needed to reproduce the cost. The script below builds a synthetic three-sub-swath layout and synthetic offset rasters at the frame-scale grid and times a single generate_insar_mask() call; the crc of the output mask lets you check current-vs-vectorized equality directly. Measured on the host described under Measurement:

implementation	wall	per pixel	mask crc
current (develop @ 0a8df45dd)	383.3 s	5.30 µs/px	0x295ba19b
vectorized (branch below)	6.7 s	0.09 µs/px (57x)	0x295ba19b

repro_insar_mask_timing.py (also kept at isce3-benchmark)

#!/usr/bin/env python3
"""Synthetic-data timing reproducer for the generate_insar_mask loop.

Requires only isce3 (with the nisar package), NumPy, GDAL and h5py — no
NISAR granules. Builds a synthetic SubSwaths layout and synthetic
geometric-coregistration offset rasters at a production-like grid size,
then times a single ``generate_insar_mask()`` call. Run it against a
pristine isce3 install to see the per-pixel-loop cost, and against the
``perf/vectorize-insar-mask`` branch to see the vectorized cost; the
printed ``mask crc`` is identical when the outputs match.

Usage:
    python3 repro_insar_mask_timing.py [lines] [samples]

Defaults to 6840 x 10581 — the RIFG interferogram grid of a NISAR
L-band frame (72.4 Mpx). Pass smaller dims for a quicker look.
"""
import sys
import tempfile
import time
import zlib

# pyre's journal reads sys.argv eagerly on import; stash it
_argv = sys.argv
sys.argv = sys.argv[:1]

import numpy as np
import h5py
from osgeo import gdal

import isce3
from nisar.products.insar.utils import generate_insar_mask

sys.argv = _argv

lines = int(sys.argv[1]) if len(sys.argv) > 1 else 6840
samples = int(sys.argv[2]) if len(sys.argv) > 2 else 10581

rng = np.random.default_rng(0)

# Three sub-swaths with jittered per-line [start, end) valid-sample
# intervals spanning the swath
edges = np.linspace(0, samples, 4).astype(np.int64)
arrays = []
for s in range(3):
    start = edges[s] + rng.integers(-20, 20, lines)
    end = edges[s + 1] + rng.integers(-20, 20, lines)
    arrays.append(np.stack([np.clip(start, 0, samples),
                            np.clip(end, 0, samples)],
                           axis=1).astype(np.int32))
subswaths = isce3.product.SubSwaths(lines, samples, arrays)


class Swath:
    def __init__(self):
        self.lines, self.samples = lines, samples

    def sub_swaths(self):
        return subswaths


class SLC:
    SwathPath = "/science/LSAR/RSLC/swaths"

    def getSwathMetadata(self, freq):
        return Swath()


tmpdir = tempfile.mkdtemp()
paths = []
drv = gdal.GetDriverByName("ENVI")
for name in ("range.off", "azimuth.off"):
    path = f"{tmpdir}/{name}"
    ds = drv.Create(path, samples, lines, 1, gdal.GDT_Float64)
    ds.GetRasterBand(1).WriteArray(rng.normal(0.0, 3.0, (lines, samples)))
    ds.FlushCache()
    ds = None
    paths.append(path)

# Empty in-memory HDF5: the inputDataExceptionMask datasets are absent,
# which exercises the zeros fallback (the per-pixel packing still runs)
h5 = h5py.File("repro", "w", driver="core", backing_store=False)

azi_idx = np.round(np.arange(lines, dtype=np.float64))
rg_idx = np.round(np.arange(samples, dtype=np.float64))

t0 = time.perf_counter()
mask = generate_insar_mask(SLC(), SLC(), h5, h5,
                           paths[0], paths[1], "A", azi_idx, rg_idx)
dt = time.perf_counter() - t0

n_px = lines * samples
print(f"generate_insar_mask: {lines}x{samples} = {n_px / 1e6:.1f} Mpx "
      f"in {dt:.1f} s ({dt / n_px * 1e6:.2f} us/px), "
      f"mask crc = {zlib.crc32(np.ascontiguousarray(mask)):#010x}")

(For full bitwise coverage — rounding edges, empty/missing sub-swath layouts, non-zero exception-mask bits — see the fixture-level verification script linked under Compatibility.)

Measurement (production-scale NISAR run)

This was found while profiling a NISAR L-band GUNW pair; the end-to-end numbers below are from that dataset. Current vs vectorized: the only treatment variable was the contents of utils.py (overlaid on an otherwise identical install tree; the run logs echo the resolved module path). One run per variant (n=1), control first then treatment, no cache reset between runs. Full logs, per-stage tables and methodology: A/B evidence.

context	current	vectorized	speedup
standalone python -m nisar.workflows.prepare_insar_hdf5 (GUNW runconfig, freq A HH)	435.4 s	86.8 s	5.02x
full insar.py GPU run, first prepare_insar_hdf5 occurrence	464.2 s	87.5 s	5.30x

In the control run, prepare_insar_hdf5 is the second-largest measured CPU-parity stage (~12% of wall, after SNAPHU unwrapping); this loop accounts for ~80% of the stage. Supporting result: in the same end-to-end GPU pair the INSAR total went 3865.0 → 3687.2 s (−177.7 s). The two prepare_insar_hdf5 occurrences together got faster by 391.1 s (−376.7 s on the first, −14.4 s on the second, ionosphere-chain one); the difference to the end-to-end delta (213.3 s) was largely offset by observed differences in unrelated stages — rubbersheet +127.7 s (a stage with documented high run-to-run variance on this host) and geo2rdr +59.9 s — with the remaining ~26 s spread across smaller per-stage differences and inter-stage glue; every other stage ratio was ≈1.0x. Since the stage is CPU/GPU parity, a similar absolute reduction is expected — not yet measured — on CPU-only runs (workflow total there ~6200 s).

Test data: NISAR L1 RSLC pair, ascending track 139 / frame 019 (Boso, Japan), 2026-07-05 / 2026-07-17, L-band DHDH, freq A HH, GUNW product type. Host: Intel Core Ultra 9 285H (16 threads), RTX 5080, NVMe scratch; isce3 v0.25.16 built from source (the function is identical on develop @ 0a8df45dd); Python 3.12, NumPy 1.26.4.

Candidate fix

Branch perf/vectorize-insar-mask (single commit 6d1fbe35b, utils.py only, +123/−52). No new C++/pybind surface is needed: the per-sub-swath per-line [start, end) interval arrays are fetched once via the existing SubSwaths.get_valid_samples_array API, and membership becomes vectorized interval tests that preserve getSampleSubSwath semantics exactly: out-of-bounds → 0, first-match-wins ordering, empty-array short-circuit, and no-sub-swath-information → 1.

Compatibility

• Bitwise-identical output (NumPy 1.26.4):
  • Fixture level: a verification script compares the vectorized function against a frozen verbatim copy of the current scalar loop on 7 fixture families covering the rounding edges (exact k + 0.5 values, negative offsets), empty and missing sub-swath layouts, out-of-swath rows/columns, out-of-bounds secondary indices, and non-zero exception-mask bits — all pass.
  • Product level: two full insar.py runs (current vs vectorized) were compared across all 558 datasets of the RIFG/RUNW/GUNW outputs; every science dataset was bitwise-identical. The only differences were processingDateTime and the repr-address runConfigurationContents, which were also observed between repeated unmodified runs.
• The current code uses two different rounding rules: int(x + 0.5) (truncation toward zero) for the sub-swath lookup and Python round() (round-half-even) for the exception-mask lookup. The vectorized code intentionally reproduces both (np.trunc(x + 0.5) / np.rint) so the output is unchanged; unifying the two rules would be a behavior change and is out of scope here.
• The exception-mask bytes are widened to uint32 before the << 16 / << 8 shifts. Besides being required for correct vectorized packing, this also addresses the NumPy ≥ 2.0 (NEP 50) hazard reported in generate_insar_mask(): inputDataExceptionMask bits silently dropped from InSAR mask layer under NumPy ≥ 2.0 (uint8 scalar left-shift overflows to 0) #335 — under NumPy 1.26.4 the old and new packing coincide, which is what the bitwise checks above pin. Those checks were produced under NumPy 1.26.4; closing generate_insar_mask(): inputDataExceptionMask bits silently dropped from InSAR mask layer under NumPy ≥ 2.0 (uint8 scalar left-shift overflows to 0) #335 itself would additionally need a direct-value packing test under NumPy ≥ 2.0, which has not been run yet.

Questions

1. Would you be open to a PR based on the linked implementation?
2. If so, would you prefer the uint32 packing change described in generate_insar_mask(): inputDataExceptionMask bits silently dropped from InSAR mask layer under NumPy ≥ 2.0 (uint8 scalar left-shift overflows to 0) #335 to be handled in the same PR?

Either way, I would include focused regression tests under a new tests/python/packages/nisar/products/insar/utils.py (registered in tests/python/packages/CMakeLists.txt), covering (a) sub-swath membership semantics (half-open intervals, first-match-wins ordering, empty and missing sub-swath layouts, out-of-bounds indices) and (b) generate_insar_mask cases pinning both rounding rules and non-zero exception-mask packing with directly asserted expected uint32 values.

---

Disclosure: this investigation and the patch were developed with assistance from AI coding tools (Claude, Codex, and Gemini). All measurements were executed on real hardware, and the evidence linked above (profiles, logs, and the bitwise verification) was generated and reviewed by the author.