r.proj: size parallel bands from a precomputed footprint grid - #7807
Draft
krcoder123 wants to merge 32 commits into
Draft
r.proj: size parallel bands from a precomputed footprint grid#7807krcoder123 wants to merge 32 commits into
krcoder123 wants to merge 32 commits into
Conversation
…re Apple M-series
Replace the whole-map RAM buffer (Path A) with a two-level band loop modeled on r.neighbors, adapted for r.proj's CRS-dependent input access: each output band's input footprint is found by back-projecting the band's edges (dense edge walk), so the loaded input strip is sized per band rather than by a fixed neighborhood stencil. Band height adapts to the memory option; if a single output row's footprint exceeds the cap (oblique or large-halo transforms) the module bails, since that case needs the tile cache path, which is not implemented here. Per-thread PROJ contexts (one PJ clone per thread) are retained. Input strips are loaded serially per band because a single fd read path is not thread-safe; each band's output is written in order. Bit-exact against the serial output at 1/2/4/8 threads, for both column-varying and row-varying inputs, with multiple bands exercised. On the test case (105.3M output cells, EPSG:4326 to EPSG:3857, nearest, memory=50) peak RSS was 130 MB versus 763 MB for the whole-map buffer. Developed with assistance from Claude (Anthropic).
The banding timers call omp_get_wtime(), which is undefined when GRASS is built without OpenMP and breaks the link in the minimum-config build. Wrap omp.h and omp_get_wtime() behind _OPENMP via a small rproj_wtime() helper that returns 0.0 without OpenMP. Also demote the PHASE_TIMERS line from G_message to G_debug, since it is benchmark scaffolding, not user output.
PROJ transformation objects are not safe for concurrent use, so a parallel module needs a private clone per thread. Add GPJ_clone_transform() and GPJ_free_transform_clone(), which bundle a cloned transform with its private PROJ context in struct gpj_transform_clone so ownership is a single unit. r.proj's parallel banding created the per-thread context with proj_context_create(), proj_clone(), and proj_destroy() directly; switch it to these helpers so the PROJ calls live in lib/proj and the module makes none.
Each band's input strip is read in parallel: read_nprocs fresh per-thread fds (Rast_open_old), a static block split of the strip rows across threads, each thread reading its disjoint rows through its own fd into its own strip slice. Rast_disable_omp_on_mask gates the parallelism (serial when a mask is present or without OpenMP); fdi remains the serial-fallback path. Env-switch choreography, compute region, write loop, band sizing, and PJ context cloning are unchanged. Experimental, not for merge.
The memory-bounded banding path halves the band height until a full-width input strip fits the memory cap. On oblique projections a single output row can back-project to an input footprint larger than the cap at any height, and that path bailed out. Add column tiling as a second search dimension. Phase 1 is identical to the current banding: halve the band height while the input strip spans the full input width, and use that result whenever a full-width band fits. Phase 2 runs only when a single full-width row still exceeds the cap. It keeps the band as tall as its output buffer allows and halves the tile width instead, so each column tile back-projects to a smaller input row span and the per-band parallel region stays populated rather than collapsing to a single row. The width search runs in two tiers to stay cheap. The upper tier estimates the worst tile strip by probing a bounded, evenly spaced subset of tiles, which is a lower bound on the true worst, and narrows to a candidate width. The lower tier validates that width with the exact per-tile edge walk and narrows further if the estimate was optimistic, so the accepted width is always exact-sized against the cap. Input strips stay full width because the raster API reads whole rows, so a tile strip is its input row span times the full input width, and tiling shrinks the row span rather than the width. Each tile loads one strip whose row span comes from the exact edge walk, one tile at a time, bounding peak memory to the worst tile rather than the whole band. Every output cell is computed once and written in row order, so the result is bit-exact with the serial output. Retain the existing fatal error only for the degenerate tile whose footprint cannot fit the cap at any width, at minimum band height; that footprint needs the tile-cache path, which is not implemented.
The tile sizing search finds, for each output band, the tallest band height and its column tiling whose input strip fits the memory cap. Neighboring bands almost always end up with the same size, since the projection changes gradually from one band to the next. The search now tries the previous band's accepted height and width first instead of restarting the descending scan from the top every time. The previous band's size is checked with the same exact edge walk acceptance test the full search uses, and the next taller height is checked to make sure it does not fit. Together these two checks confirm the reused size is the tallest fitting answer, the same result the full scan would have returned. If either check fails, the code falls back to the full descending search, so the worst case costs the same as before. Because acceptance is decided by the same test in both paths, the resulting bands and tiles are identical to before and the output is bit for bit unchanged. In the common case a band is sized in two edge walks instead of a full descending scan.
The memory-bounded band sizing walks each output tile's perimeter to bound the range of input rows the tile needs, then loads that strip. A tile whose interior contains a geographic pole has its northmost or southmost latitude at the pole, in the tile interior, where the perimeter walk never samples it. The strip was therefore sized too small, and projecting a pole-containing map aborted with a "Band strip under-sized" error at every memory setting, though the projection itself was well defined. This computes, once per map, each geographic pole that lies within the input's latitude coverage: its coordinate in the output projection and its input row. When an output tile's rectangle contains a pole, that pole's input row is folded into the tile's row span, so the height and width search sees the true footprint and shrinks pole tiles until they fit the memory cap. The loaded strip then covers every row the fill reads. Only lat/lon input is handled, where a pole is at latitude 90 or -90. If the pole's coordinate transform fails or returns a non-finite value the pole is skipped and the existing under-size guard stays as the backstop. A map with no pole in the output frame is unaffected: the row spans, the band and tile partition, and the output are byte for byte unchanged.
The banded compute path dispatched a nearest-only strip reader for every resampling method, so bilinear, bicubic, lanczos and their fallback variants silently produced nearest-neighbor output. Add interp_strip.c with strip counterparts of the cache kernels (strip_bilinear, strip_cubic, strip_lanczos, and the three _f fallbacks). They read the in-RAM full-width FCELL band strip the banded path already loads, using the same base index, bounds, weights, and null fallback as the readcell-cache kernels in bilinear.c, cubic.c, and lanczos.c. A strip_kernels[] table, ordered like menu[], resolves each method to its strip counterpart once after option parsing; nearest keeps the existing interpolate_strip reader in slot 0. Output is bitwise identical to serial r.proj for all seven methods across the test datasets.
Serial r.proj computes each output row's northing by subtracting ns_res row by row. The banded path computed it directly as north - ns_res/2 - row * ns_res, which can differ by one ulp when ns_res is not exactly representable. Nearest is unaffected, but for the other methods the shifted interpolation weights changed a few cells (29 of 76M on the EPSG:3035 test). Precompute the row northings once with the serial recurrence and use that array in both the sizing walk and the compute loop.
The banded path aborted when the memory cap could not hold even one output row's input strip. In practice that needs an input wider than about cap/20 columns; poles and oblique projections do not trigger it. Warn with the minimum memory that keeps the parallel path, then finish the run through the old serial readcell cache. Output is bitwise identical to serial r.proj. Failed transforms set the cell NULL like the banded path does. R_PROJ_FORCE_TILECACHE forces the fallback for testing.
GPJ_clone_transform did not check the results of proj_context_create() and proj_clone(). Either can return NULL on failure, and the NULL would otherwise surface later as a crash deep inside PROJ when the cloned transform is first used. Both are now checked and fail with G_fatal_error naming the call. r.proj calls this once per worker thread, so a failure terminates the process from inside the parallel region; that is intended, since a clone failure leaves the thread with no usable transform.
…input A polar-stereographic output frame centered on a pole, reprojecting an input truncated below that pole (e.g. input reaching 89 degrees, not 90), aborted with "Band strip under-sized" (or, before that guard existed, silently read garbage). The banded strip sizing sampled only the output tile's perimeter, so the frame-center-proximal interior cell that reaches the input's northernmost edge row was never seen, and the strip loaded too few input rows. The pole footprint fold already handled a pole lying inside the input map; it now also folds in the input's edge row (0 or rows-1) when a pole outside the input's latitude coverage still projects into the frame. Gated by the existing point-in-rect test, so bands that do not image a pole compute byte-identical spans; verified unchanged on non-pole frames.
Adds pytest tests that verify r.proj's banded parallel output matches its serial output on generated CI-sized data. r.proj has no nprocs option, so each run sets OMP_NUM_THREADS on a per-call environment copy to select the serial (1) or parallel (N) path without mutating shared state. Four tests: bilinear identity (with a nearest-vs-bilinear dispatch-liveness guard so a silent fallback to nearest cannot pass the check vacuously), nearest identity under a constrained memory cap that forces band sizing, nearest identity into a pole-centered frame, and a forced tile-cache fallback (R_PROJ_FORCE_TILECACHE) compared against the banded path to cover both code paths. Inputs are integer CELL below 2^24 so the FCELL readcell cache round-trips losslessly.
Ports the 7 interpolation-method tests from testsuite/test_rproj.py to a parametrized pytest (r_proj_methods_test.py) over generated CI-sized data, with references captured from the serial binary. Removes the gunittest file. Drops its 4 output-format tests (test_list_output_plain/json, test_print_output_plain/json), which asserted NC-SPM-specific -l/-p output; r.proj's -l and -p flags are left without test coverage as a result.
Overlapping input rows are kept between single-tile bands and only new rows are read, instead of re-reading each band's full span. Cuts the input read phase about 70 percent on the wide LAEA benchmark; output is bitwise identical to serial.
The band output buffer is now double-buffered when running with more than one thread: one thread writes the previous band's rows in order while the rest compute the current band. Single-thread runs keep the sequential write. Hides most of the output write time at higher thread counts; output is bitwise identical to serial.
The method reference tests and the gunittest-to-pytest testsuite replacement move to a dedicated PR (branch fix-rproj-tests) so this PR stays focused on the OpenMP parallelization and that PR catches serial regressions on its own. This drops raster/r.proj/tests/r_proj_methods_test.py and restores raster/r.proj/testsuite/test_rproj.py to its main state; the split PR carries the deletion, so this restore is temporary and goes away once that PR merges. The parallel-correctness tests (r_proj_parallel_test.py) and their conftest.py stay here; conftest.py gains a note that it is duplicated on the split branch.
Add the standard G_OPT_M_NPROCS option so the compute thread count can be set with nprocs= instead of only OMP_NUM_THREADS. A value above zero overrides OMP_NUM_THREADS and zero keeps the OpenMP default. The option is read once through compute_nprocs() before the band fit search, so it drives the compute region, the per-thread read fds, and the output double-buffer together. The parallel-correctness tests now pass nprocs= instead of setting OMP_NUM_THREADS. Also shorten the main.c comments to flowing prose, dropping restated design narration and internal shorthand while keeping the load-bearing rationale.
Add a benchmark that sweeps the nprocs= thread count from 1 to 8 at two memory caps and plots the time, speedup, and efficiency metrics, following the r.param.scale benchmark template with grass.benchmark. It builds a source project and reprojects a generated raster from EPSG:4326 into EPSG:3857 in a temporary database, so it is self-contained.
Add a footprint grid of input row spans, one grid row per output row and thirty-two column blocks wide, in a boundary-sampled and a column-exact variant. The grid only observes in this commit. The live fit search still steers every band height and tile width, and the grid span is compared against the search span for each rectangle the search evaluates. The comparison is gated by the R_PROJ_FG_VERIFY environment variable, so a normal run builds no grid and prints nothing. Each boundary cell carries a one row sampling margin, since the column samples can miss a curve between them by a fraction of a row. The margin is applied after the two variants are compared, so the variant report still measures the raw sampling error.
Band heights now come from a scan of the footprint grid instead of the Phase-1 halving search, which is removed along with its seed machinery. Phase 2 and the rest of the pipeline are unchanged, and the output stays bitwise identical to the serial reference.
Tile widths now come from a scan of the footprint grid, with tiles built as runs of whole column blocks and a single over-cap check that routes to the serial fallback. The former tile-width search is removed and the output stays bitwise identical to the serial reference.
The parallel nearest path now rounds each input value through a 32-bit float as the strip is read, the same step the serial cache-based read uses. Outputs then match the serial reference for every input, including integers above 2^24 where the float step changes the value.
The footprint grid was validated against the perimeter-walk search while it was being brought up. That verification path, its environment flag, and the counters it printed are no longer needed, and the runtime under-size check remains as the backstop.
The nearest strip reader is renamed to strip_nearest to match its sibling kernels, the footprint grid pointer to band_grid now that there is only one grid, and the pole row field to pole_row.
Shorten the longer comments to plain descriptions of what the code does and drop the internal jargon. No code changes.
The band loop co-sizes each band's height and tile width from the footprint grid. It grows the height while a full-width band fits the memory cap and turns to whole column tiles only when even one full-width row does not fit, so an output row whose columns project across a wide span of input rows is read in tall bands instead of one row at a time. The output values are identical to the serial result.
github-actions
Bot
added
raster
Shortens the docstrings and comments in the tests and the benchmark script.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This draft PR is to show a second way to parallelize r.proj. It’s still pretty similar to #7627 but the main difference is #7627 sized bands by prediction. It guessed a band height, tested whether the input fit the memory cap, halved the guess and tried again. It also had a second search of the same kind for tile widths. This was around 340 lines of guessing in total.
This PR replaces the sizing with a table. When the run starts, the module projects the output edges once and writes down where everything lands in the input. For example, Row 0, block 0 lands on these input rows, it writes that down. Row 0, block 1 lands here, it writes that down, so on and so forth.
Once the table is filled, sizing is just reading the table. A band starts at full width rows and grows one at a time until the memory cap can’t take anymore rows. If even one full width row cannot fit, the band is read as column tiles instead. The tile starts wide and drops one block at a time until it fits. So the module never predicts. It just reads the table, stops when the cap is reached, and moves to the next band. When nothing fits at all it knows up front and falls back to the serial tile cache.
This is around 58 lines smaller, with about 340 lines of sizing search replaced by the grid. I also fixed up a bunch of the comments. These are the benchmarks for the various scenarios I ran:
Benchmarks for all tested scenarios
"easy" EPSG:4326 to EPSG:3857, memory=50, method=nearest
Serial: 10.26s (measured same session)
"easy" EPSG:4326 to EPSG:3857, memory=300, method=nearest
Serial: 9.56s (measured same session)
"easy" EPSG:4326 to EPSG:3857, memory=50, method=bilinear
Serial: 10.34s (measured same session)
"easy" EPSG:4326 to EPSG:3857, memory=300, method=bilinear
Serial: 10.18s (measured same session)
"laea" EPSG:4326 to EPSG:3035, memory=50, method=nearest
Serial: 12.85s (measured same session)
"laea" EPSG:4326 to EPSG:3035, memory=300, method=nearest
Serial: 10.95s (measured same session)
"laea" EPSG:4326 to EPSG:3035, memory=50, method=bilinear
Serial: 13.40s (measured same session)
"laea" EPSG:4326 to EPSG:3035, memory=300, method=bilinear
Serial: 11.07s (measured same session)
"pole" EPSG:4326 to EPSG:3413, memory=50, method=nearest
Serial: 30.17s (measured same session)
"pole" EPSG:4326 to EPSG:3413, memory=300, method=nearest
Serial: 29.30s (measured same session)
"pole" EPSG:4326 to EPSG:3413, memory=50, method=bilinear
Serial: 30.37s (measured same session)
"pole" EPSG:4326 to EPSG:3413, memory=300, method=bilinear
Serial: 29.48s (measured same session)
"easy" EPSG:4326 to EPSG:3857, memory=50, method=lanczos
Serial: 16.53s (measured same session)
"laea" EPSG:4326 to EPSG:3035, memory=50, method=lanczos
Serial: 18.22s (measured same session)
Fallback sweep: memory-cap bail boundary (out_mult=2, N=8)