Skip to content

Vectorize the per-pixel loop in generate_insar_mask (#354) - #359

Open
s-sasaki-earthsea-wizard wants to merge 2 commits into
isce-framework:developfrom
s-sasaki-earthsea-wizard:perf/vectorize-insar-mask
Open

Vectorize the per-pixel loop in generate_insar_mask (#354)#359
s-sasaki-earthsea-wizard wants to merge 2 commits into
isce-framework:developfrom
s-sasaki-earthsea-wizard:perf/vectorize-insar-mask

Conversation

@s-sasaki-earthsea-wizard

Copy link
Copy Markdown

Summary

Implements the change proposed in #354: generate_insar_mask() built the InSAR mask layer with a pure-Python double loop over every output pixel (two scalar pybind SubSwaths.get_sample_sub_swath calls, scalar indexing and a list.append per pixel). The per-pixel work is elementwise with no loop-carried dependency, so this PR replaces the inner loop with NumPy array operations built on the existing SubSwaths.get_valid_samples_array bulk API. On the NISAR L-band frame-scale grid (72.4 Mpx) the function drops from 383.3 s to 6.7 s (57x) on synthetic data, and the containing prepare_insar_hdf5 stage runs 5x faster with bitwise-identical output on a real production run (435.4 s to 86.8 s). Full measurement evidence — py-spy profile attribution, the synthetic self-contained reproducer, and the production A/B — is in #354.

Fixes #335, closes #354.

Changes

  • generate_insar_mask() (python/packages/nisar/products/insar/utils.py): fetch each sub-swath's per-line [start, end) valid-sample interval array once via SubSwaths.get_valid_samples_array and evaluate sub-swath membership as vectorized interval tests in a new _subswath_numbers helper, preserving the scalar get_sample_sub_swath semantics: out-of-bounds -> 0, first-match-wins ordering, empty-array short-circuit, and no-sub-swath-information -> 1. Both rounding rules are preserved exactly: int(x + 0.5) (truncation toward zero) for the sub-swath lookup via np.trunc, and Python round() (round-half-even) for the exception-mask lookup via np.rint. The GDAL row reads are unchanged.
  • The inputDataExceptionMask bytes are widened to uint32 before the << 16 / << 8 packing shifts. Under NumPy >= 2 (NEP 50) the previous uint8-scalar shifts overflow to 0 and silently drop those bits (generate_insar_mask(): inputDataExceptionMask bits silently dropped from InSAR mask layer under NumPy ≥ 2.0 (uint8 scalar left-shift overflows to 0) #335); the explicit widening produces the intended packing under both promotion regimes.
  • New regression tests (tests/python/packages/nisar/products/insar/utils.py, registered in the test CMakeLists): _subswath_numbers is checked against the scalar SubSwaths.get_sample_sub_swath API as oracle, and generate_insar_mask against a per-pixel scalar reference on synthetic fixtures covering adversarial offsets (exact k + 0.5 half-integers where the two rounding rules diverge, large out-of-swath pushes), differing reference/secondary dimensions, empty and absent sub-swath layouts, the missing-inputDataExceptionMask path, exact uint32 packing of exception bytes with the MSB set (the generate_insar_mask(): inputDataExceptionMask bits silently dropped from InSAR mask layer under NumPy ≥ 2.0 (uint8 scalar left-shift overflows to 0) #335 regression case), and a negative secondary index where truncation toward zero and floor diverge.

Behavior

Under NumPy 1.x the output is bitwise-identical to the current implementation — verified on synthetic fixtures covering the semantic edges (standalone harness linked from #354) and on a production run where all mask layers across RIFG/RUNW/GUNW match bitwise. The one intended difference: under NumPy >= 2, exception-mask bytes are packed correctly instead of being silently dropped, which is the #335 fix.

Testing

  • The new test file passes (11 tests) in a Python 3.12 / NumPy 1.26 / GDAL 3.12 environment against this branch.
  • The seven behavior tests (differential, bit packing, rounding) also pass unchanged against the current develop scalar implementation under NumPy 1.26 — they encode the pre-existing semantics, not the new implementation's, and both implementations satisfy them.
  • Synthetic timing reproducer and production-scale A/B: see Performance: vectorizing generate_insar_mask() makes prepare_insar_hdf5 5x faster #354.

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.

generate_insar_mask built the InSAR mask layer with a pure-Python
double loop over every output pixel, calling the scalar
SubSwaths.get_sample_sub_swath twice per pixel. On a NISAR L-band
frame (RIFG interferogram grid 6840x10581) this costs ~4.2 us/px,
~300 s per product, and makes prepare_insar_hdf5 the largest
non-GPU-addressable stage of the InSAR workflow (~560 s, CPU/GPU
parity since the cost is interpreter-bound).

Replace the inner loop with numpy array operations:

- Fetch each sub-swath's per-line [start, end) valid-sample interval
  array once via the existing SubSwaths.get_valid_samples_array API
  and evaluate membership as vectorized interval tests
  (_subswath_numbers), preserving get_sample_sub_swath semantics:
  out-of-bounds -> 0, first-match-wins ordering, empty-array
  short-circuit, and no-sub-swath-information -> 1.
- Preserve both rounding rules exactly: int(x + 0.5) (truncation
  toward zero) for the sub-swath lookup via np.trunc, and Python
  round() (round-half-even) for the exception-mask lookup via
  np.rint.
- Widen the exception-mask bytes to uint32 before the << 16 / << 8
  shifts. Under NumPy >= 2.0 (NEP 50) the previous uint8-scalar
  shifts overflow to 0 and silently drop those bits (isce-framework#335); the
  explicit widening produces the intended packing under both
  promotion regimes.

The GDAL row reads and the produced mask values are unchanged; output
is bitwise-identical to the previous implementation (verified against
a frozen copy of the scalar loop on fixtures covering the rounding
edges, empty/missing sub-swath layouts, out-of-swath rows/columns,
out-of-bounds secondary indices, and high exception-mask bits, under
NumPy 1.26).
Cover the semantics the vectorized implementation must preserve:

- _subswath_numbers against the scalar
  SubSwaths.get_sample_sub_swath API as oracle: out-of-bounds -> 0,
  first-match-wins ordering, the empty valid-samples-array
  short-circuit, and the no-sub-swath-information -> 1 path;
- generate_insar_mask against a per-pixel scalar reference
  (_compute_subswath_mask_id plus Python-int bit packing) on synthetic
  fixtures with adversarial offsets (exact k + 0.5 half-integers,
  large out-of-swath pushes), differing reference/secondary
  dimensions, empty and absent sub-swath layouts, and the missing
  inputDataExceptionMask dataset path;
- exact uint32 packing of exception-mask bytes with the MSB set,
  which the previous uint8-scalar << 16 / << 8 shifts silently drop
  under NumPy >= 2 scalar promotion (isce-framework#335);
- the rounding asymmetry between the sub-swath lookup (int(x + 0.5),
  truncation toward zero) and the exception-mask lookup
  (round-half-even), including a negative secondary index where
  truncation toward zero and floor diverge.

The seven behavior tests pass unchanged against the pre-vectorization
scalar implementation under NumPy 1.26, confirming they encode the
existing semantics rather than the new implementation's.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

1 participant