diff --git a/README.md b/README.md index 3275579..6859417 100644 --- a/README.md +++ b/README.md @@ -62,8 +62,8 @@ RasterMap( backend=backend, band="state", color_map={0: "#ffffff", 1: "#2f8f6e"}, - labels={0: "morta", 1: "viva"}, - title=f"Padrões clássicos sobre fronteiras de bloco", + labels={0: "dead", 1: "alive"}, + title="Classic patterns over block boundaries", ) env.run() @@ -191,8 +191,8 @@ is what tends to expose it. ## Two input paths: GeoTIFF/VRT or Zarr -`geotiff_io.py`/`mosaic_io.py` (TIFF/VRT, via `rasterio`) and -`zarr_io.py` (via `zarr`) are **interchangeable** input paths — both +`disk/io/geotiff.py` (TIFF/VRT, via `rasterio`) and +`disk/io/zarr.py` (via `zarr`) are **interchangeable** input paths — both populate the same `MemmapRasterWorkspace`, block by block, without materializing the full array in RAM. Swap the loader; the rest of the pipeline (halo, disk, models) doesn't change. @@ -250,7 +250,7 @@ never imported by `haloexec`'s runtime code). ## Loading GeoTIFF straight to disk -`geotiff_io.py` (`load_geotiff_into_workspace`) loads a real GeoTIFF +`disk/io/geotiff.py` (`load_geotiff_into_workspace`) loads a real GeoTIFF block by block directly into a `MemmapRasterWorkspace`, via `rasterio.windows.Window` — never materializing a whole band in RAM. It uses the `band_spec` convention already established by @@ -281,11 +281,10 @@ routing, watershed delineation — where a cell's value can, in principle, depend on the entire domain, not just its immediate neighbors. -Generalized from a domain-specific student prototype (tidal -connectivity via `scipy.ndimage.binary_propagation`), though that -prototype's specific rule is not part of this primitive — only the -orchestration pattern was extracted: a small halo plus repeated global -sweeps until no block changes, instead of one large halo. Each sweep +The primitive holds only the orchestration pattern — the rule itself +(e.g. tidal connectivity via `scipy.ndimage.binary_propagation`) is +supplied by the caller: a small halo plus repeated global sweeps until +no block changes, instead of one large halo. Each sweep writes its result **immediately** back (Gauss-Seidel, via `write_block_core_in_place` — no ping-pong), so a block processed later in the same sweep already sees the update from a block processed @@ -346,16 +345,13 @@ seeds). ## Disk layer (grids larger than RAM) -`disk_backend.py` (`MemmapRasterWorkspace`) and `disk_sync_model.py` +`disk/workspace.py` (`MemmapRasterWorkspace`) and `disk/sync_model.py` (`DiskChunkedSyncRasterModel`) generalize the same domain decomposition for grids that don't fit in memory, using `np.memmap` with double-buffering and checkpointing. **No dependency on dissmodel** in -`disk_backend.py` — reusable by any framework. +`disk/workspace.py` — reusable by any framework. -Extracted and generalized from a student prototype (a domain-specific -preprocessing pipeline) that already had correct memmap+double-buffer -+checkpoint mechanics, but tied to domain-specific state names. Here -arrays are named generically (a `name -> dtype` dict), with no coupling +Arrays are named generically (a `name -> dtype` dict), with no coupling to any particular domain's variable names. ```python @@ -425,7 +421,7 @@ The cost is low — a typical workspace has ~1500 blocks, so the index is on the order of 1500 bits, versus the ~371 MB a per-cell mask would cost. `read_block_*` for an absent block would return the declared `nodata` without touching disk. The point of attention is -`disk_sync_model.py`, which today writes every block unconditionally +`disk/sync_model.py`, which today writes every block unconditionally in `write_block_core` — the change would need to decide whether a fully-nodata block should be written back at all. @@ -438,7 +434,7 @@ evidence that the axis order stored on disk isn't guaranteed. Reproduced with real `xarray`, writing exactly the way disscube's `VariableWriter` writes (`da.to_dataset(...).to_zarr(...)`): a **square** array with `(x, y)` axes instead of `(y, x)` has the -**same shape** in both cases — the shape check in `zarr_io.py` didn't +**same shape** in both cases — the shape check in `disk/io/zarr.py` didn't catch the inversion. Without a fix, this corrupted row/column **silently**, with no error at all. diff --git a/docs/README.md b/docs/README.md index 519958c..405a12d 100644 --- a/docs/README.md +++ b/docs/README.md @@ -125,10 +125,10 @@ from haloexec import ( DiskChunkedRasterCellularAutomaton, ) -# DiskChunkedRasterCellularAutomaton é um mixin cooperativo: precisa vir -# primeiro no MRO, com uma classe real do dissmodel (aqui, -# RasterCellularAutomaton) como segunda base — do contrário setup() -# não tem para onde delegar via super() e a classe falha ao instanciar. +# DiskChunkedRasterCellularAutomaton is a cooperative mixin: it must come +# first in the MRO, with a real dissmodel class (here, +# RasterCellularAutomaton) as the second base — otherwise setup() has +# nowhere to delegate through super() and the class fails to instantiate. class LargeScaleGameOfLife(DiskChunkedRasterCellularAutomaton, RasterCellularAutomaton): def rule(self, arrays: dict[str, np.ndarray]) -> dict[str, np.ndarray]: state = arrays["state"] diff --git a/docs/api_reference.md b/docs/api_reference.md index aee3092..a55e158 100644 --- a/docs/api_reference.md +++ b/docs/api_reference.md @@ -41,7 +41,7 @@ This document provides a comprehensive, exhaustive reference for all public clas ## `haloexec.engine` -Primitivas de decomposição de domínio sem dependências externas. +Domain decomposition primitives with no external dependencies. ### `Block` ```python @@ -94,7 +94,7 @@ Resolves the external ghost cell fill value for a given array. ## `haloexec.disk.workspace` -Gerenciador de arrays bidimensionais em disco baseados em `np.memmap` com suporte a double-buffering e checkpoints. +Manager of two-dimensional on-disk arrays based on `np.memmap`, with double-buffering and checkpoints. ### `HaloWindow` ```python diff --git a/docs/tutorials_and_recipes.md b/docs/tutorials_and_recipes.md index 157969d..4e306e5 100644 --- a/docs/tutorials_and_recipes.md +++ b/docs/tutorials_and_recipes.md @@ -72,9 +72,9 @@ from dissmodel.core import Environment from dissmodel.geo.raster.cellular_automaton import RasterCellularAutomaton from haloexec import MemmapRasterWorkspace, DiskChunkedRasterCellularAutomaton -# Mixin cooperativo: precisa de RasterCellularAutomaton como segunda -# base (MRO) — herança simples aqui falha em runtime (setup() não tem -# para onde delegar via super()). +# Cooperative mixin: needs RasterCellularAutomaton as the second base +# (MRO) — single inheritance here fails at runtime (setup() has nowhere +# to delegate through super()). class MassiveGameOfLife(DiskChunkedRasterCellularAutomaton, RasterCellularAutomaton): def rule(self, arrays: dict[str, np.ndarray]) -> dict[str, np.ndarray]: state = arrays["state"] diff --git a/examples/gol/gol_patterns_haloexec.py b/examples/gol/gol_patterns_haloexec.py index 96a56c5..02aafaf 100644 --- a/examples/gol/gol_patterns_haloexec.py +++ b/examples/gol/gol_patterns_haloexec.py @@ -1,23 +1,23 @@ """ -Game of Life com padrões clássicos + haloexec + RasterMap -============================================================ -Posiciona padrões conhecidos (glider, blinker, beacon, toad, block, -pulsar) DELIBERADAMENTE sobre as fronteiras de bloco, para que seja -fácil verificar visualmente se o halo está sincronizando corretamente -(um padrão que atravessa a borda de um bloco e continua se comportando -como deveria é a evidência visual mais direta). - -Requisitos ----------- +Game of Life with classic patterns + haloexec + RasterMap +========================================================= +Places well-known patterns (glider, blinker, beacon, toad, block, +pulsar) DELIBERATELY over block boundaries, so it is easy to check +visually whether the halo synchronizes correctly (a pattern that +crosses a block edge and keeps behaving as it should is the most direct +visual evidence). + +Requirements +------------ pip install dissmodel - pip install -e /caminho/para/dissmodel-ca - pip install -e /caminho/para/haloexec + pip install -e /path/to/dissmodel-ca + pip install -e /path/to/haloexec -Uso ---- - python gol_patterns_haloexec.py +Usage +----- + python examples/gol/gol_patterns_haloexec.py -Sem display interativo, os PNGs caem em ./raster_map_frames/. +Without an interactive display, the PNGs go to ./raster_map_frames/. """ from __future__ import annotations @@ -31,7 +31,7 @@ # --------------------------------------------------------------------------- -# Modelo: GameOfLife com halo +# Model: GameOfLife with a halo # --------------------------------------------------------------------------- class GameOfLifeHalo(HaloChunkedRasterCellularAutomaton): def setup( @@ -59,7 +59,7 @@ def rule(self, arrays: dict) -> dict: def place(grid: np.ndarray, pattern: list[list[int]], top: int, left: int) -> None: - """Escreve um padrão na grade a partir de (top, left).""" + """Write a pattern into the grid starting at (top, left).""" arr = np.array(pattern) h, w = arr.shape grid[top:top + h, left:left + w] = arr @@ -75,19 +75,19 @@ def place(grid: np.ndarray, pattern: list[list[int]], top: int, left: int) -> No grid = np.zeros((ROWS, COLS), dtype=np.int8) -# Posicionados DE PROPÓSITO sobre fronteiras de bloco (múltiplos de 10) -# — o pior caso para testar a sincronização do halo. -# CORRIGIDO: beacon movido de col 25 -> col 33 (na versão original, -# beacon em (9,25) 4x4 colidia com pulsar em (4,20) 13x13 -- a área -# do beacon [9:13, 25:29] cai inteira dentro da área do pulsar -# [4:17, 20:33]. place() sobrescreve por atribuição direta, então os -# dois padrões ficariam corrompidos silenciosamente, sem erro nenhum. -place(grid, PATTERNS["glider"], 8, 8) # atravessa o cruzamento (10,10) na diagonal -place(grid, PATTERNS["blinker"], 20, 5) # atravessa a borda horizontal em r=20 -place(grid, PATTERNS["beacon"], 9, 33) # atravessa a borda horizontal em r=10 -place(grid, PATTERNS["toad"], 29, 15) # atravessa a borda horizontal em r=30 -place(grid, PATTERNS["block"], 19, 19) # atravessa o cruzamento de 4 blocos (10,20)+(10,20) -place(grid, PATTERNS["pulsar"], 4, 20) # oscilador maior, período 3 +# Placed ON PURPOSE over block boundaries (multiples of 10) +# — the worst case for testing halo synchronization. +# FIXED: beacon moved from col 25 -> col 33 (in the original version, +# the 4x4 beacon at (9,25) collided with the 13x13 pulsar at (4,20) -- +# the beacon's area [9:13, 25:29] falls entirely inside the pulsar's +# area [4:17, 20:33]. place() overwrites by direct assignment, so both +# patterns would have been silently corrupted, with no error at all. +place(grid, PATTERNS["glider"], 8, 8) # crosses the (10,10) corner diagonally +place(grid, PATTERNS["blinker"], 20, 5) # crosses the horizontal edge at r=20 +place(grid, PATTERNS["beacon"], 9, 33) # crosses the horizontal edge at r=10 +place(grid, PATTERNS["toad"], 29, 15) # crosses the horizontal edge at r=30 +place(grid, PATTERNS["block"], 19, 19) # sits on the corner where 4 blocks meet (20,20) +place(grid, PATTERNS["pulsar"], 4, 20) # larger oscillator, period 3 backend = raster_grid(rows=ROWS, cols=COLS, attrs={"state": grid}) @@ -101,14 +101,14 @@ def place(grid: np.ndarray, pattern: list[list[int]], top: int, left: int) -> No ) # --------------------------------------------------------------------------- -# Visualização +# Visualization # --------------------------------------------------------------------------- RasterMap( backend=backend, band="state", color_map={0: "#ffffff", 1: "#2f8f6e"}, - labels={0: "morta", 1: "viva"}, - title=f"Padrões clássicos sobre fronteiras de bloco (halo={HALO})", + labels={0: "dead", 1: "alive"}, + title=f"Classic patterns over block boundaries (halo={HALO})", ) # --------------------------------------------------------------------------- diff --git a/examples/gol/scale_30m.py b/examples/gol/scale_30m.py new file mode 100644 index 0000000..97c503d --- /dev/null +++ b/examples/gol/scale_30m.py @@ -0,0 +1,151 @@ +""" +Scale test: ~30 million pixels, writing the TIFF in windows (never +materializing the whole grid in RAM), loading it block by block, and +running dissmodel_ca's real GameOfLife through +DiskChunkedRasterCellularAutomaton (haloexec) -- the same scale +challenge a real coastal model faces. + +Compares against a monolithic reference AT THE SAME SCALE (30M uint8 +cells fit in RAM as a single array, ~30 MB -- what must not happen is +materializing while WRITING/LOADING the large file, which is what this +script actually tests). +""" + +import time +from pathlib import Path + +import numpy as np +import rasterio +from dissmodel.core import Environment +from dissmodel.geo import raster_grid +from dissmodel_ca.models.game_of_life_raster import GameOfLife +from rasterio.transform import from_origin +from rasterio.windows import Window + +from haloexec import DiskChunkedRasterCellularAutomaton, MemmapRasterWorkspace, load_geotiff_into_workspace + + +class GameOfLifeHalo(DiskChunkedRasterCellularAutomaton, GameOfLife): + pass + + +def memory_mb() -> dict: + """RssAnon (real heap) from /proc/self/status -- the metric that + proves materialization, not VmRSS/ru_maxrss (which include the page + cache of mapped files, always high with memmap without meaning a + problem).""" + values = {} + with open("/proc/self/status") as f: + for line in f: + for key in ("VmRSS", "RssAnon", "RssFile"): + if line.startswith(key + ":"): + values[key] = int(line.split()[1]) / 1024 + return values + + +def write_tiff_in_windows(path: Path, height: int, width: int, density: float, + seed: int, block: int = 512) -> None: + """Write the GeoTIFF block by block through rasterio.windows.Window + -- never allocates the whole (height, width) grid in RAM at once. + Deterministic RNG per block position (reproducible).""" + transform = from_origin(500_000.0, 9_700_000.0, 30.0, 30.0) + master_rng = np.random.default_rng(seed) + + with rasterio.open( + str(path), "w", driver="GTiff", height=height, width=width, + count=1, dtype="uint8", crs="EPSG:31984", transform=transform, + tiled=True, blockxsize=block, blockysize=block, compress="lzw", + ) as dst: + for r0 in range(0, height, block): + r1 = min(r0 + block, height) + for c0 in range(0, width, block): + c1 = min(c0 + block, width) + block_seed = int(master_rng.integers(0, 2**31 - 1)) ^ (r0 * 92821 + c0) + rng = np.random.default_rng(block_seed & 0xFFFFFFFF) + data = (rng.random((r1 - r0, c1 - c0)) < density).astype("uint8") + dst.write(data, 1, window=Window(c0, r0, c1 - c0, r1 - r0)) + + +def main(): + # ~30 million pixels + HEIGHT, WIDTH = 5480, 5480 # 30,030,400 cells + DENSITY = 0.35 + SEED = 42 + GENERATIONS = 5 + BLOCK = 256 + HALO = 1 + + tmp = Path("/tmp/haloexec_scale_30m") + tmp.mkdir(exist_ok=True) + tif_path = tmp / "initial_state.tif" + + print(f"Grid: {HEIGHT}x{WIDTH} = {HEIGHT*WIDTH:,} cells " + f"(~{HEIGHT*WIDTH/1024**2:.1f} MB per uint8 array)") + + # ── 1. write the TIFF in windows ──────────────────────────────── + rss_before = memory_mb() + t0 = time.time() + write_tiff_in_windows(tif_path, HEIGHT, WIDTH, DENSITY, SEED, block=BLOCK) + t_write = time.time() - t0 + rss_after_write = memory_mb() + print(f"\n[1/3] TIFF written in {t_write:.1f}s -- " + f"RssAnon={rss_after_write['RssAnon']:.1f}MB " + f"(delta since start: {rss_after_write['RssAnon']-rss_before['RssAnon']:.1f}MB)") + + # ── 2. load block by block into the workspace ─────────────────── + t0 = time.time() + ws = MemmapRasterWorkspace.create( + root=tmp / "workspace", shape=(HEIGHT, WIDTH), + arrays={"state": np.uint8}, block_h=BLOCK, block_w=BLOCK, halo=HALO, + ) + load_geotiff_into_workspace(ws, tif_path, [("state", "uint8", 0)]) + t_load = time.time() - t0 + rss_after_load = memory_mb() + print(f"[2/3] Loaded in {t_load:.1f}s -- " + f"RssAnon={rss_after_load['RssAnon']:.1f}MB " + f"(delta since writing: {rss_after_load['RssAnon']-rss_after_write['RssAnon']:.1f}MB)") + + # ── 3. run GameOfLife on disk+halo, through dissmodel ─────────── + t0 = time.time() + env = Environment(start_time=1, end_time=GENERATIONS) + GameOfLifeHalo(workspace=ws, halo=HALO, boundary_value=0) + env.run() + ws.flush() + t_run = time.time() - t0 + rss_after_run = memory_mb() + print(f"[3/3] {GENERATIONS} generations in {t_run:.1f}s " + f"({t_run/GENERATIONS*1000:.0f}ms/generation) -- " + f"RssAnon={rss_after_run['RssAnon']:.1f}MB " + f"(delta since loading: {rss_after_run['RssAnon']-rss_after_load['RssAnon']:.1f}MB)") + + disk_result = ws.snapshot("state") + + print("\n=== memory summary ===") + print(f"size of one full array: {HEIGHT*WIDTH/1024**2:.1f} MB") + print(f"final RssAnon: {rss_after_run['RssAnon']:.1f} MB " + f"(ratio to 1 array: {rss_after_run['RssAnon']/(HEIGHT*WIDTH/1024**2):.2f}x)") + print(f"final RssFile: {rss_after_run['RssFile']:.1f} MB (page cache, not materialization)") + + # ── 4. equivalence against a monolithic reference AT THE SAME SCALE ── + # 30M uint8 cells fit in RAM as a SINGLE array (~30 MB) -- what must + # not happen is materializing while WRITING and LOADING the file, + # which was shown above through RssAnon. + print("\n[extra] building a monolithic reference at the same scale for the equivalence proof...") + with rasterio.open(str(tif_path)) as ds: + state0 = ds.read(1) # here we DO materialize, on purpose, only for the golden reference + backend_mono = raster_grid(rows=HEIGHT, cols=WIDTH, attrs={"state": state0.copy()}) + env_mono = Environment(start_time=1, end_time=GENERATIONS) + GameOfLife(backend=backend_mono) + t0 = time.time() + env_mono.run() + t_mono = time.time() - t0 + golden = backend_mono.arrays["state"].copy() + + n_diff = int(np.sum(golden != disk_result)) + print(f"monolithic run took {t_mono:.1f}s ({t_mono/GENERATIONS*1000:.0f}ms/generation)") + print(f"\ndisk-vs-monolithic differences: {n_diff} of {HEIGHT*WIDTH:,} cells") + print(f"IDENTICAL? {n_diff == 0}") + + +if __name__ == "__main__": + main() diff --git a/examples/gol/teste_escala_com_rastermap.py b/examples/gol/scale_30m_rastermap.py similarity index 54% rename from examples/gol/teste_escala_com_rastermap.py rename to examples/gol/scale_30m_rastermap.py index eedc44b..d4cfe75 100644 --- a/examples/gol/teste_escala_com_rastermap.py +++ b/examples/gol/scale_30m_rastermap.py @@ -1,10 +1,10 @@ """ -Teste de escala com visualização controlada via RasterMap: -30 milhões de pixels em disco + halo, salvando quadros PNG apenas em -gerações/anos selecionados para não sobrecarregar a memória nem o tempo de CPU. +Scale test with controlled visualization through RasterMap: +30 million pixels on disk + halo, saving PNG frames only at selected +generations/years so neither memory nor CPU time is overloaded. -Uso: - python examples/gol_patterns/teste_escala_com_rastermap.py +Usage: + python examples/gol/scale_30m_rastermap.py """ from __future__ import annotations @@ -33,21 +33,21 @@ class GameOfLifeHalo(DiskChunkedRasterCellularAutomaton, GameOfLife): pass -def memoria_mb() -> dict[str, float]: - """RssAnon (heap real) via /proc/self/status.""" - valores = {} +def memory_mb() -> dict[str, float]: + """RssAnon (real heap) from /proc/self/status.""" + values = {} with open("/proc/self/status") as f: for line in f: - for chave in ("VmRSS", "RssAnon", "RssFile"): - if line.startswith(chave + ":"): - valores[chave] = int(line.split()[1]) / 1024 - return valores + for key in ("VmRSS", "RssAnon", "RssFile"): + if line.startswith(key + ":"): + values[key] = int(line.split()[1]) / 1024 + return values -def gerar_tiff_em_janelas( +def write_tiff_in_windows( path: Path, height: int, width: int, density: float, seed: int, block: int = 512 ) -> None: - """Gera GeoTIFF escrevendo bloco a bloco via Window (nunca aloca grade inteira em RAM).""" + """Write a GeoTIFF block by block through Window (never allocates the whole grid in RAM).""" transform = from_origin(500_000.0, 9_700_000.0, 30.0, 30.0) master_rng = np.random.default_rng(seed) @@ -62,35 +62,35 @@ def gerar_tiff_em_janelas( c1 = min(c0 + block, width) block_seed = int(master_rng.integers(0, 2**31 - 1)) ^ (r0 * 92821 + c0) rng = np.random.default_rng(block_seed & 0xFFFFFFFF) - dado = (rng.random((r1 - r0, c1 - c0)) < density).astype("uint8") - dst.write(dado, 1, window=Window(c0, r0, c1 - c0, r1 - r0)) + data = (rng.random((r1 - r0, c1 - c0)) < density).astype("uint8") + dst.write(data, 1, window=Window(c0, r0, c1 - c0, r1 - r0)) def main() -> None: - HEIGHT, WIDTH = 5480, 5480 # ~30.030.400 células + HEIGHT, WIDTH = 5480, 5480 # ~30,030,400 cells DENSITY = 0.35 SEED = 42 GENERATIONS = 5 BLOCK = 256 HALO = 1 - # Anos / passos a salvar na visualização - ANOS_PARA_SALVAR = [1, 3, 5] + # Years / steps to save in the visualization + STEPS_TO_SAVE = [1, 3, 5] - tmp = Path("/tmp/teste_escala_rastermap") + tmp = Path("/tmp/haloexec_scale_rastermap") shutil.rmtree(tmp, ignore_errors=True) tmp.mkdir(parents=True, exist_ok=True) - tif_path = tmp / "estado_inicial.tif" + tif_path = tmp / "initial_state.tif" - print("=== Teste de Escala com Visualização por Checkpoints ===") - print(f"Grade: {HEIGHT}x{WIDTH} = {HEIGHT*WIDTH:,} células (~{HEIGHT*WIDTH/1024**2:.1f} MB por array uint8)") - print(f"Salvando quadros PNG apenas nos passos: {ANOS_PARA_SALVAR}\n") + print("=== Scale test with checkpoint visualization ===") + print(f"Grid: {HEIGHT}x{WIDTH} = {HEIGHT*WIDTH:,} cells (~{HEIGHT*WIDTH/1024**2:.1f} MB per uint8 array)") + print(f"Saving PNG frames only at steps: {STEPS_TO_SAVE}\n") - # 1. Gerar TIFF em janelas + # 1. Write the TIFF in windows t0 = time.time() - gerar_tiff_em_janelas(tif_path, HEIGHT, WIDTH, DENSITY, SEED, block=BLOCK) - print(f"[1/4] TIFF inicial gerado em janelas ({time.time() - t0:.1f}s)") + write_tiff_in_windows(tif_path, HEIGHT, WIDTH, DENSITY, SEED, block=BLOCK) + print(f"[1/4] Initial TIFF written in windows ({time.time() - t0:.1f}s)") - # 2. Criar workspace e carregar + # 2. Create the workspace and load t0 = time.time() ws = MemmapRasterWorkspace.create( root=tmp / "workspace", @@ -101,41 +101,41 @@ def main() -> None: halo=HALO, ) load_geotiff_into_workspace(ws, tif_path, [("state", "uint8", 0)]) - print(f"[2/4] Carregado no workspace em disco ({time.time() - t0:.1f}s)") + print(f"[2/4] Loaded into the on-disk workspace ({time.time() - t0:.1f}s)") - # 3. Configurar adapter e visualizador - # stride=4 decima a grade de 5480x5480 para 1370x1370 para plotagem super rápida + # 3. Set up the adapter and the viewer + # stride=4 decimates the 5480x5480 grid to 1370x1370 for fast plotting backend_adapter = WorkspaceRasterBackend(ws, stride=4) env = Environment(start_time=1, end_time=GENERATIONS) GameOfLifeHalo(workspace=ws, halo=HALO, boundary_value=0) - # CheckpointRasterMap só executa nos passos indicados + # CheckpointRasterMap only runs at the given steps CheckpointRasterMap( backend=backend_adapter, band="state", color_map={0: "#ffffff", 1: "#2f8f6e"}, - labels={0: "morta", 1: "viva"}, - title="Game of Life 30M (Disco+Halo)", + labels={0: "dead", 1: "alive"}, + title="Game of Life 30M (Disk+Halo)", save_frames=True, - save_steps=ANOS_PARA_SALVAR, + save_steps=STEPS_TO_SAVE, auto_mask=False, ) - # 4. Rodar simulação - rss_antes = memoria_mb() + # 4. Run the simulation + rss_before = memory_mb() t0 = time.time() env.run() ws.flush() t_exec = time.time() - t0 - rss_depois = memoria_mb() + rss_after = memory_mb() - print(f"\n[3/4] Simulação finalizada em {t_exec:.1f}s ({t_exec/GENERATIONS*1000:.0f}ms/geração)") - print(f"RssAnon final: {rss_depois['RssAnon']:.1f} MB (delta: {rss_depois['RssAnon'] - rss_antes['RssAnon']:.1f} MB)") + print(f"\n[3/4] Simulation finished in {t_exec:.1f}s ({t_exec/GENERATIONS*1000:.0f}ms/generation)") + print(f"final RssAnon: {rss_after['RssAnon']:.1f} MB (delta: {rss_after['RssAnon'] - rss_before['RssAnon']:.1f} MB)") out_dir = Path("raster_map_frames") pngs = sorted(out_dir.glob("state_step_*.png")) if out_dir.exists() else [] - print(f"\n[4/4] Quadros gerados em {out_dir}/:") + print(f"\n[4/4] Frames written to {out_dir}/:") for p in pngs: print(f" - {p.name} ({p.stat().st_size / 1024:.1f} KB)") diff --git a/examples/gol/teste_escala_30m.py b/examples/gol/teste_escala_30m.py deleted file mode 100644 index 7c3df83..0000000 --- a/examples/gol/teste_escala_30m.py +++ /dev/null @@ -1,150 +0,0 @@ -""" -Teste de escala: ~30 milhões de pixels, gerando o TIFF em janelas -(nunca materializa a grade inteira em RAM), carregando bloco a bloco, -e rodando o GameOfLife real do dissmodel_ca via -DiskChunkedRasterCellularAutomaton (haloexec) -- o mesmo desafio de -escala que o modelo real do manguezal vai enfrentar. - -Compara contra uma referência monolítica NA MESMA ESCALA (30M células -uint8 cabem em RAM como um único array, ~30MB -- o que não cabe é o -processo de GERAR/CARREGAR via arquivo grande sem materializar, que é -o que este script realmente testa). -""" - -import time -from pathlib import Path - -import numpy as np -import rasterio -from dissmodel.core import Environment -from dissmodel.geo import raster_grid -from dissmodel_ca.models.game_of_life_raster import GameOfLife -from rasterio.transform import from_origin -from rasterio.windows import Window - -from haloexec import DiskChunkedRasterCellularAutomaton, MemmapRasterWorkspace, load_geotiff_into_workspace - - -class GameOfLifeHalo(DiskChunkedRasterCellularAutomaton, GameOfLife): - pass - - -def memoria_mb() -> dict: - """RssAnon (heap real) via /proc/self/status -- a métrica que - prova materialização, não VmRSS/ru_maxrss (inclui cache de página - de arquivo mapeado, sempre alto em memmap sem indicar problema).""" - valores = {} - with open("/proc/self/status") as f: - for line in f: - for chave in ("VmRSS", "RssAnon", "RssFile"): - if line.startswith(chave + ":"): - valores[chave] = int(line.split()[1]) / 1024 - return valores - - -def gerar_tiff_em_janelas(path: Path, height: int, width: int, density: float, - seed: int, block: int = 512) -> None: - """Gera o GeoTIFF escrevendo bloco a bloco via rasterio.windows.Window - -- nunca aloca a grade (height, width) inteira em RAM de uma vez. - RNG determinística por posição de bloco (reprodutível).""" - transform = from_origin(500_000.0, 9_700_000.0, 30.0, 30.0) - master_rng = np.random.default_rng(seed) - - with rasterio.open( - str(path), "w", driver="GTiff", height=height, width=width, - count=1, dtype="uint8", crs="EPSG:31984", transform=transform, - tiled=True, blockxsize=block, blockysize=block, compress="lzw", - ) as dst: - for r0 in range(0, height, block): - r1 = min(r0 + block, height) - for c0 in range(0, width, block): - c1 = min(c0 + block, width) - block_seed = int(master_rng.integers(0, 2**31 - 1)) ^ (r0 * 92821 + c0) - rng = np.random.default_rng(block_seed & 0xFFFFFFFF) - dado = (rng.random((r1 - r0, c1 - c0)) < density).astype("uint8") - dst.write(dado, 1, window=Window(c0, r0, c1 - c0, r1 - r0)) - - -def main(): - # ~30 milhões de pixels - HEIGHT, WIDTH = 5480, 5480 # 30.030.400 células - DENSITY = 0.35 - SEED = 42 - GENERATIONS = 5 - BLOCK = 256 - HALO = 1 - - tmp = Path("/tmp/teste_escala_30m") - tmp.mkdir(exist_ok=True) - tif_path = tmp / "estado_inicial.tif" - - print(f"Grade: {HEIGHT}x{WIDTH} = {HEIGHT*WIDTH:,} células " - f"(~{HEIGHT*WIDTH/1024**2:.1f} MB por array uint8)") - - # ── 1. gerar o TIFF em janelas ────────────────────────────────── - rss_antes = memoria_mb() - t0 = time.time() - gerar_tiff_em_janelas(tif_path, HEIGHT, WIDTH, DENSITY, SEED, block=BLOCK) - t_geracao = time.time() - t0 - rss_pos_geracao = memoria_mb() - print(f"\n[1/3] TIFF gerado em {t_geracao:.1f}s -- " - f"RssAnon={rss_pos_geracao['RssAnon']:.1f}MB " - f"(delta desde inicio: {rss_pos_geracao['RssAnon']-rss_antes['RssAnon']:.1f}MB)") - - # ── 2. carregar bloco a bloco no workspace ────────────────────── - t0 = time.time() - ws = MemmapRasterWorkspace.create( - root=tmp / "workspace", shape=(HEIGHT, WIDTH), - arrays={"state": np.uint8}, block_h=BLOCK, block_w=BLOCK, halo=HALO, - ) - load_geotiff_into_workspace(ws, tif_path, [("state", "uint8", 0)]) - t_carga = time.time() - t0 - rss_pos_carga = memoria_mb() - print(f"[2/3] Carregado em {t_carga:.1f}s -- " - f"RssAnon={rss_pos_carga['RssAnon']:.1f}MB " - f"(delta desde geracao: {rss_pos_carga['RssAnon']-rss_pos_geracao['RssAnon']:.1f}MB)") - - # ── 3. rodar GameOfLife em disco+halo, via dissmodel ──────────── - t0 = time.time() - env = Environment(start_time=1, end_time=GENERATIONS) - GameOfLifeHalo(workspace=ws, halo=HALO, boundary_value=0) - env.run() - ws.flush() - t_execucao = time.time() - t0 - rss_pos_execucao = memoria_mb() - print(f"[3/3] {GENERATIONS} gerações em {t_execucao:.1f}s " - f"({t_execucao/GENERATIONS*1000:.0f}ms/geração) -- " - f"RssAnon={rss_pos_execucao['RssAnon']:.1f}MB " - f"(delta desde carga: {rss_pos_execucao['RssAnon']-rss_pos_carga['RssAnon']:.1f}MB)") - - resultado_disco = ws.snapshot("state") - - print("\n=== resumo de memoria ===") - print(f"tamanho de um array completo: {HEIGHT*WIDTH/1024**2:.1f} MB") - print(f"RssAnon final: {rss_pos_execucao['RssAnon']:.1f} MB " - f"(razao sobre 1 array: {rss_pos_execucao['RssAnon']/(HEIGHT*WIDTH/1024**2):.2f}x)") - print(f"RssFile final: {rss_pos_execucao['RssFile']:.1f} MB (cache de pagina, nao e materializacao)") - - # ── 4. equivalencia contra referencia monolitica NA MESMA ESCALA ── - # 30M celulas uint8 cabem em RAM como array UNICO (~30MB) -- o que - # nao cabe/nao deveria ser feito e materializar durante GERACAO e - # CARGA do arquivo, que ja foi provado acima via RssAnon. - print("\n[extra] gerando referencia monolitica na mesma escala para prova de equivalencia...") - with rasterio.open(str(tif_path)) as ds: - estado0 = ds.read(1) # aqui SIM materializamos, de proposito, so para a referencia golden - backend_mono = raster_grid(rows=HEIGHT, cols=WIDTH, attrs={"state": estado0.copy()}) - env_mono = Environment(start_time=1, end_time=GENERATIONS) - GameOfLife(backend=backend_mono) - t0 = time.time() - env_mono.run() - t_mono = time.time() - t0 - golden = backend_mono.arrays["state"].copy() - - n_diff = int(np.sum(golden != resultado_disco)) - print(f"monolitico rodou em {t_mono:.1f}s ({t_mono/GENERATIONS*1000:.0f}ms/geracao)") - print(f"\ndivergencias disco-vs-monolitico: {n_diff} de {HEIGHT*WIDTH:,} celulas") - print(f"IDENTICO? {n_diff == 0}") - - -if __name__ == "__main__": - main() diff --git a/pyproject.toml b/pyproject.toml index 1ac2896..2f10996 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ build-backend = "setuptools.build_meta" [project] name = "haloexec" version = "0.1.0" -description = "Motor genérico de execução de Autômatos Celulares por decomposição de domínio (blocos+halo)" +description = "Generic cellular automaton execution engine using domain decomposition (blocks + halo)" readme = "README.md" license = "MIT" license-files = ["LICENSE"] @@ -38,43 +38,42 @@ docs = [ "mkdocs-material>=9.5", "markdown-callouts>=0.4", # github-callouts extension (> [!WARNING]) ] -# Opcional: só necessário para HaloChunkedRasterCellularAutomaton, -# HaloChunkedSyncRasterModel, DiskChunkedSyncRasterModel -# (dissmodel_ca.py, sync_model.py, disk_sync_model.py) -- os únicos 3 -# módulos do haloexec que dependem de dissmodel. Os outros 5 -# (engine.py, disk_backend.py, convergence.py, geotiff_io.py, -# zarr_io.py) funcionam sem dissmodel instalado. +# Optional: only needed for the dissmodel adapters +# (HaloChunkedRasterCellularAutomaton, HaloChunkedSyncRasterModel, +# DiskChunkedRasterCellularAutomaton, DiskChunkedSyncRasterModel, in +# ram/ and disk/*_model.py / disk/cellular_automaton.py). engine.py, +# disk/workspace.py, disk/convergence.py and disk/io/ work without it. dissmodel = ["dissmodel>=0.6.3"] -# Opcional: só necessário para geotiff_io.py (carregamento de GeoTIFF/VRT -# direto para MemmapRasterWorkspace, bloco a bloco). +# Optional: only needed for disk/io/geotiff.py (loading GeoTIFF/VRT +# straight into MemmapRasterWorkspace, block by block). geotiff = ["rasterio>=1.3"] -# Opcional: só necessário para zarr_io.py (carregamento de Zarr direto -# para MemmapRasterWorkspace — segunda opção de entrada, pensada para -# consumir DerivedVariable do disscube). Requer zarr-python 3 (API -# create_array/dimension_names), que por sua vez exige Python >= 3.11. +# Optional: only needed for disk/io/zarr.py (loading Zarr straight into +# MemmapRasterWorkspace — the second input path, meant for disscube's +# DerivedVariable). Requires zarr-python 3 (create_array/dimension_names +# API), which in turn requires Python >= 3.11. zarr = ["zarr>=3"] -# Opcional: só necessário para tests/test_zarr_axis_order_regression.py, -# que reproduz o padrão real de escrita do disscube (VariableWriter usa -# xarray.to_zarr()) para validar a correção de ordem de eixos. +# Optional: only needed for tests/test_zarr_axis_order_regression.py, +# which reproduces disscube's real write pattern (VariableWriter uses +# xarray.to_zarr()) to validate the axis-order fix. zarr-test = ["zarr>=3", "xarray>=2024.10"] -# Opcional: só necessário para tests/test_geomosaic_integration.py, que -# prova que load_geotiff_into_workspace lê corretamente um mosaico -# produzido pelo geomosaic (pacote separado, sem relação de runtime -# com haloexec — só usado aqui como prova de integração em teste). +# Optional: only needed for tests/test_geomosaic_integration.py, which +# proves that load_geotiff_into_workspace correctly reads a mosaic +# produced by geomosaic (a separate package with no runtime relation to +# haloexec — used here only as an integration proof in tests). geomosaic = ["geomosaic @ git+https://github.com/LambdaGeo/geomosaic.git"] -# Opcional: só necessário para examples/gol_patterns/gol_patterns_haloexec.py -# e tests/test_gol_patterns_example.py (o teste de equivalência desse -# exemplo) -- dissmodel-ca fornece PATTERNS (padrões clássicos de Game -# of Life) e dissmodel.visualization.RasterMap. +# Optional: only needed for examples/gol/ and +# tests/test_gol_patterns_example.py (the equivalence test for that +# example) -- dissmodel-ca provides PATTERNS (classic Game of Life +# patterns) and GameOfLife. dissmodel-ca = ["dissmodel-ca @ git+https://github.com/DisSModel/dissmodel-ca.git"] [tool.setuptools.packages.find] where = ["src"] [tool.pytest.ini_options] -# Restringe a descoberta padrão do pytest a tests/ -- sem isso, os -# scripts de exemplo em examples/ seriam coletados automaticamente por -# um `pytest` sem argumento na raiz do repositório. +# Restrict pytest's default discovery to tests/ -- otherwise the example +# scripts in examples/ would be collected by a bare `pytest` at the +# repository root. testpaths = ["tests"] [tool.ruff] diff --git a/scripts/generate_and_benchmark.py b/scripts/generate_and_benchmark.py index fc253cc..958d6ea 100644 --- a/scripts/generate_and_benchmark.py +++ b/scripts/generate_and_benchmark.py @@ -1,24 +1,24 @@ """ -Gera uma grade sintética grande DIRETO EM DISCO (bloco a bloco, nunca -materializada inteira em RAM) e roda Game of Life via -MemmapRasterWorkspace, medindo tempo de execução e pico de RSS de RAM. +Generate a large synthetic grid STRAIGHT ON DISK (block by block, never +materialized whole in RAM) and run Game of Life through +MemmapRasterWorkspace, measuring run time and peak RAM RSS. -Objetivo: provar empiricamente que o footprint de memória fica -limitado ao tamanho do bloco (+halo), não ao tamanho da grade — -independentemente de quão grande for a grade em disco. +Goal: show empirically that the memory footprint stays bounded by the +block size (+halo), not by the grid size — however large the grid on +disk is. -Uso ---- +Usage +----- python scripts/generate_and_benchmark.py --shape 20000 20000 \\ --block 512 512 --halo 1 --generations 5 --density 0.35 \\ --root /tmp/haloexec_bench python scripts/generate_and_benchmark.py --shape 5000 5000 \\ - --block 128 128 --generations 20 --root /tmp/bench_pequeno + --block 128 128 --generations 20 --root /tmp/bench_small -Troque a regra (_game_of_life_rule) por qualquer outra função -`dict[str, np.ndarray] -> dict[str, np.ndarray]` para testar outro -modelo simples — a mecânica de geração/benchmark não muda. +Swap the rule (_game_of_life_rule) for any other +`dict[str, np.ndarray] -> dict[str, np.ndarray]` function to test +another simple model — the generation/benchmark mechanics stay the same. """ from __future__ import annotations @@ -49,8 +49,8 @@ def _game_of_life_rule(padded: dict[str, np.ndarray], halo: int = 1) -> dict[str def generate_synthetic_on_disk( ws: MemmapRasterWorkspace, name: str, density: float, seed: int ) -> None: - """Popula um array do workspace bloco a bloco, com RNG determinística - por bloco — nunca aloca a grade inteira em RAM de uma vez.""" + """Fill one workspace array block by block, with a deterministic RNG + per block — never allocates the whole grid in RAM at once.""" master_rng = np.random.default_rng(seed) for block in ws.blocks(): block_seed = int(master_rng.integers(0, 2**31 - 1)) ^ (block.r0 * 92821 + block.c0) @@ -61,14 +61,14 @@ def generate_synthetic_on_disk( def memory_breakdown_mb() -> dict[str, float]: - """Quebra de RSS via /proc/self/status: RssAnon é o que o processo - de fato alocou no heap (arrays Python/numpy mantidos vivos); - RssFile é cache de páginas de arquivos mapeados (mmap) tocadas — - reclamável pelo kernel sob pressão de memória, NÃO é o mesmo que - "a grade inteira está materializada no processo". ru_maxrss/VmRSS - somam os dois, o que é enganoso para workflows baseados em mmap: - RssFile cresce com o volume de dados TOCADO ao longo do tempo - (cumulativo), não com o quanto está retido de uma vez.""" + """RSS breakdown from /proc/self/status: RssAnon is what the process + actually allocated on the heap (Python/numpy arrays kept alive); + RssFile is the page cache of touched memory-mapped files — + reclaimable by the kernel under memory pressure, NOT the same as + "the whole grid is materialized in the process". ru_maxrss/VmRSS + add the two, which is misleading for mmap-based workflows: RssFile + grows with the volume of data TOUCHED over time (cumulative), not + with how much is held at once.""" values = {} with open("/proc/self/status") as f: for line in f: @@ -93,11 +93,11 @@ def run_benchmark( if root.exists(): shutil.rmtree(root) - grid_bytes = shape[0] * shape[1] # uint8: 1 byte/célula - print(f"Grade: {shape[0]}x{shape[1]} = {shape[0]*shape[1]:,} células " - f"(~{grid_bytes / 1024**2:.1f} MB por array, x2 slots x2 discos " - f"= ~{grid_bytes * 4 / 1024**2:.1f} MB em disco)") - print(f"Bloco: {block_h}x{block_w}, halo={halo}, gerações={generations}") + grid_bytes = shape[0] * shape[1] # uint8: 1 byte/cell + print(f"Grid: {shape[0]}x{shape[1]} = {shape[0]*shape[1]:,} cells " + f"(~{grid_bytes / 1024**2:.1f} MB per array, x2 slots x2 disks " + f"= ~{grid_bytes * 4 / 1024**2:.1f} MB on disk)") + print(f"Block: {block_h}x{block_w}, halo={halo}, generations={generations}") rss_before = memory_breakdown_mb() @@ -107,9 +107,9 @@ def run_benchmark( block_h=block_h, block_w=block_w, halo=halo, ) generate_synthetic_on_disk(ws, "state", density=density, seed=seed) - t_geracao = time.time() - t0 + t_write = time.time() - t0 m = memory_breakdown_mb() - print(f"Geração sintética em disco: {t_geracao:.2f}s | " + print(f"Synthetic generation on disk: {t_write:.2f}s | " f"RssAnon={m['RssAnon']:.1f}MB RssFile={m['RssFile']:.1f}MB " f"VmRSS={m['VmRSS']:.1f}MB") @@ -122,40 +122,40 @@ def run_benchmark( ws.swap_buffers() ws.checkpoint(step) ws.flush() - t_execucao = time.time() - t0 - - m_depois = memory_breakdown_mb() - print(f"Execução ({generations} gerações): {t_execucao:.2f}s " - f"({t_execucao/generations*1000:.1f} ms/geração)") - print(f"RssAnon (heap real do processo): {m_depois['RssAnon']:.1f} MB " - f"(delta desde o início: {m_depois['RssAnon'] - rss_before['RssAnon']:.1f} MB)") - print(f"RssFile (cache de páginas mmap, reclamável): {m_depois['RssFile']:.1f} MB") - print(f"VmRSS total (soma dos dois, é o que ru_maxrss mediria): {m_depois['VmRSS']:.1f} MB") - print(f"Tamanho de um array completo em disco: {grid_bytes / 1024**2:.1f} MB") - print("→ RssAnon é a métrica correta para 'quanto o processo materializou " - "de fato'; RssFile cresce com o volume TOCADO acumulado (cache), " - "não com o que está retido de uma vez.") + t_run = time.time() - t0 + + m_after = memory_breakdown_mb() + print(f"Run ({generations} generations): {t_run:.2f}s " + f"({t_run/generations*1000:.1f} ms/generation)") + print(f"RssAnon (the process's real heap): {m_after['RssAnon']:.1f} MB " + f"(delta since start: {m_after['RssAnon'] - rss_before['RssAnon']:.1f} MB)") + print(f"RssFile (mmap page cache, reclaimable): {m_after['RssFile']:.1f} MB") + print(f"total VmRSS (sum of the two, what ru_maxrss would measure): {m_after['VmRSS']:.1f} MB") + print(f"Size of one full array on disk: {grid_bytes / 1024**2:.1f} MB") + print("→ RssAnon is the right metric for 'how much the process actually " + "materialized'; RssFile grows with the accumulated volume TOUCHED " + "(cache), not with what is held at once.") if not keep: shutil.rmtree(root) - print(f"Workspace removido ({root}). Use --keep para preservar.") + print(f"Workspace removed ({root}). Use --keep to keep it.") else: - print(f"Workspace preservado em {root}") + print(f"Workspace kept at {root}") if __name__ == "__main__": parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--shape", type=int, nargs=2, default=[20000, 20000], - metavar=("ALTURA", "LARGURA")) + metavar=("HEIGHT", "WIDTH")) parser.add_argument("--block", type=int, nargs=2, default=[512, 512], - metavar=("BLOCO_H", "BLOCO_W")) + metavar=("BLOCK_H", "BLOCK_W")) parser.add_argument("--halo", type=int, default=1) parser.add_argument("--generations", type=int, default=5) parser.add_argument("--density", type=float, default=0.35) parser.add_argument("--seed", type=int, default=42) parser.add_argument("--root", type=Path, default=Path("/tmp/haloexec_bench")) parser.add_argument("--keep", action="store_true", - help="não apagar o workspace ao final") + help="do not delete the workspace at the end") args = parser.parse_args() run_benchmark( diff --git a/src/haloexec/__init__.py b/src/haloexec/__init__.py index 496455e..098d375 100644 --- a/src/haloexec/__init__.py +++ b/src/haloexec/__init__.py @@ -23,12 +23,11 @@ "sweep_until_convergence", ] -# Os adaptadores dissmodel (HaloChunkedRasterCellularAutomaton, -# HaloChunkedSyncRasterModel, DiskChunkedSyncRasterModel) são opcionais -# -- os módulos acima funcionam sem dissmodel instalado. Só ficam -# disponíveis se o extra "dissmodel" estiver instalado -# (pip install "haloexec[dissmodel]"). Mesmo padrão usado em -# pymangue/__init__.py para CMMAModel. +# The dissmodel adapters (HaloChunkedRasterCellularAutomaton, +# HaloChunkedSyncRasterModel, DiskChunkedSyncRasterModel, +# DiskChunkedRasterCellularAutomaton) are optional -- the modules above +# work without dissmodel installed. They are only available when the +# "dissmodel" extra is installed (pip install "haloexec[dissmodel]"). try: from .disk.cellular_automaton import DiskChunkedRasterCellularAutomaton from .disk.sync_model import DiskChunkedSyncRasterModel, workspace_arrays_for_sync_model diff --git a/src/haloexec/disk/backend.py b/src/haloexec/disk/backend.py index ff2559a..3e6bf7b 100644 --- a/src/haloexec/disk/backend.py +++ b/src/haloexec/disk/backend.py @@ -1,9 +1,9 @@ """ -Adaptador de RasterBackend para MemmapRasterWorkspace. +RasterBackend adapter for MemmapRasterWorkspace. -Permite que componentes do ecossistema dissmodel (como visualizadores -RasterMap, exportadores ou coletores) acessem os arrays do slot de leitura -atual do workspace sem materializar nem duplicar a grade inteira em RAM. +Lets components of the dissmodel ecosystem (such as RasterMap viewers, +exporters or collectors) access the arrays of the workspace's current +read slot without materializing or duplicating the whole grid in RAM. """ from __future__ import annotations @@ -15,23 +15,24 @@ class WorkspaceRasterBackend: """ - Adaptador leve que expõe um MemmapRasterWorkspace com a interface de - RasterBackend (shape e dicionário de arrays). + Lightweight adapter that exposes a MemmapRasterWorkspace through the + RasterBackend interface (shape and a dictionary of arrays). - Os arrays expostos são referências diretas (np.memmap) ao slot de leitura - atual do workspace, respeitando double-buffering e swap_buffers. + The exposed arrays are direct references (np.memmap) to the + workspace's current read slot, honouring double-buffering and + swap_buffers. - Parâmetros + Parameters ---------- workspace : MemmapRasterWorkspace - Workspace em disco a ser adaptado. + On-disk workspace to adapt. stride : int, default=1 - Fator de subamostragem espacial (decimação) ao acessar arrays. - Útil para visualização em larga escala (ex.: 30M+ células), - reduzindo significativamente o tempo de renderização e o uso de memória - do matplotlib sem alterar os dados no disco. - nodata_value : float | int | None, default=None - Valor de nodata opcional para integração com extent masks do RasterMap. + Spatial subsampling (decimation) factor applied when accessing + arrays. Useful for large-scale visualization (e.g. 30M+ cells), + greatly reducing matplotlib's rendering time and memory use + without changing the data on disk. + nodata_value : float | None, default=None + Optional nodata value, for RasterMap's extent masks. """ def __init__( @@ -48,7 +49,7 @@ def __init__( @property def arrays(self) -> dict[str, np.ndarray]: - """Dicionário com os arrays do slot de leitura corrente.""" + """Dictionary with the arrays of the current read slot.""" slot = self.workspace.checkpoint_data["read_slot"] memmaps = self.workspace._slots[slot] if self.stride == 1: @@ -56,11 +57,11 @@ def arrays(self) -> dict[str, np.ndarray]: return {name: mm[::self.stride, ::self.stride] for name, mm in memmaps.items()} def get(self, name: str) -> np.ndarray: - """Retorna o array correspondente ao nome.""" + """Return the array with the given name.""" return self.arrays[name] def snapshot(self) -> dict[str, np.ndarray]: - """Retorna uma cópia dos arrays do slot de leitura atual.""" + """Return a copy of the arrays of the current read slot.""" return {k: np.asarray(v).copy() for k, v in self.arrays.items()} def __repr__(self) -> str: diff --git a/src/haloexec/disk/cellular_automaton.py b/src/haloexec/disk/cellular_automaton.py index d3bc647..80ea85d 100644 --- a/src/haloexec/disk/cellular_automaton.py +++ b/src/haloexec/disk/cellular_automaton.py @@ -1,20 +1,20 @@ """ -Integração com dissmodel: DiskChunkedRasterCellularAutomaton. - -Equivalente em disco de ram/cellular_automaton.py::HaloChunkedRasterCellularAutomaton -— mesmo contrato de rule() (`rule(arrays) -> dict`), mas lendo/escrevendo -via MemmapRasterWorkspace em vez de materializar a grade inteira em RAM. - -Preenche uma lacuna real: existia adaptador de disco pra modelos -SyncRasterModel (disk/sync_model.py::DiskChunkedSyncRasterModel), mas -não para modelos RasterCellularAutomaton — qualquer AC escrito com -rule() (como dissmodel_ca.models.game_of_life_raster.GameOfLife) só -rodava em RAM até este módulo existir. - -Uso: `class GameOfLifeHalo(DiskChunkedRasterCellularAutomaton, GameOfLife): pass` -reusa o rule() real da classe original sem reescrever nada — mesmo -princípio de composição cooperativa (MRO) já usado em -disk/sync_model.py com FloodModel/MangroveModel. +dissmodel integration: DiskChunkedRasterCellularAutomaton. + +Disk counterpart of ram/cellular_automaton.py::HaloChunkedRasterCellularAutomaton +— same rule() contract (`rule(arrays) -> dict`), but reading/writing +through MemmapRasterWorkspace instead of materializing the whole grid +in RAM. + +Complements disk/sync_model.py::DiskChunkedSyncRasterModel (the disk +adapter for SyncRasterModel models) for RasterCellularAutomaton models: +any CA written with rule() (such as +dissmodel_ca.models.game_of_life_raster.GameOfLife) can run from disk. + +Usage: `class GameOfLifeHalo(DiskChunkedRasterCellularAutomaton, GameOfLife): pass` +reuses the original class's real rule() without rewriting anything — +the same cooperative-composition (MRO) principle used in +disk/sync_model.py with FloodModel/MangroveModel. """ from __future__ import annotations @@ -26,11 +26,11 @@ class DiskChunkedRasterCellularAutomaton: """ - Mixin que processa rule() de um RasterCellularAutomaton em blocos - lidos de um MemmapRasterWorkspace, sem nunca materializar a grade - inteira em RAM. + Mixin that runs a RasterCellularAutomaton's rule() in blocks read + from a MemmapRasterWorkspace, never materializing the whole grid in + RAM. - Ordem de herança (MRO): este mixin deve vir primeiro, ex.: + Inheritance order (MRO): this mixin must come first, e.g. `class GameOfLifeHalo(DiskChunkedRasterCellularAutomaton, GameOfLife)`. """ @@ -46,10 +46,10 @@ def setup( self.halo = workspace.halo if halo is None else halo self.boundary_value = boundary_value - # Placeholder leve: RasterBackend(shape=...) não aloca arrays, - # só existe para satisfazer o contrato de RasterModel.setup() - # (self.backend = backend; self.shape = backend.shape). O - # backend real por bloco é criado dentro de execute(). + # Lightweight placeholder: RasterBackend(shape=...) allocates no + # arrays; it only satisfies RasterModel.setup()'s contract + # (self.backend = backend; self.shape = backend.shape). The real + # per-block backend is created inside execute(). placeholder = RasterBackend(shape=workspace.shape) super().setup(backend=placeholder, state_attr=state_attr, **kwargs) @@ -68,7 +68,7 @@ def execute(self) -> None: self.backend = block_backend self.shape = block_backend.shape try: - updates = self.rule(block_backend.snapshot()) # mesmo contrato de sempre + updates = self.rule(block_backend.snapshot()) # the usual contract finally: self.backend = real_backend self.shape = real_shape diff --git a/src/haloexec/disk/convergence.py b/src/haloexec/disk/convergence.py index 788780c..27a920d 100644 --- a/src/haloexec/disk/convergence.py +++ b/src/haloexec/disk/convergence.py @@ -1,39 +1,36 @@ """ -Varreduras repetidas até convergência — para problemas de dependência -ESPACIAL NÃO-LIMITADA (conectividade, roteamento de fluxo, delineação -de bacia), onde nenhum halo de tamanho fixo resolve sozinho: o valor -de uma célula pode depender, em princípio, do domínio inteiro. - -Generalizado de um protótipo de aluno -(chunked_engine.py::propagar_conectividade), que resolvia exatamente -esse problema para conectividade de maré via scipy.ndimage.binary_propagation -por bloco, repetindo varreduras globais até nenhum bloco mudar mais -nada. A regra específica dele (binary_propagation) NÃO faz parte desta -primitiva — só o padrão de orquestração (halo pequeno + repetição, em -vez de halo grande de uma vez) foi extraído. Qualquer `rule` que -opere sobre uma janela com halo e devolva um núcleo atualizado serve. - -Por que isso é o padrão certo para dependência não-limitada ------------------------------------------------------------------ -Um halo de tamanho fixo só deixa informação atravessar UMA fronteira -de bloco por chamada. Repetir a varredura N vezes (N = número de -blocos, no pior caso) deixa a "frente" de propagação andar um bloco -por rodada — depois de N rodadas, necessariamente alcançou qualquer -célula alcançável no domínio inteiro, não importa o tamanho. É o -mesmo princípio de BFS distribuído por rounds limitados. - -Gauss-Seidel, não Jacobi ------------------------------------------------------------------ -Cada bloco escreve o resultado IMEDIATAMENTE de volta no mesmo slot de -leitura (via write_block_core_in_place, sem ping-pong) — um bloco -processado depois, na MESMA varredura, já enxerga a atualização de um -bloco processado antes. Isso acelera a convergência (menos varreduras -necessárias) sem comprometer o resultado final, desde que a regra seja -um operador monótono (ex.: conectividade só cresce, nunca encolhe) — -nesse caso o ponto fixo final independe da ordem dos blocos, só o -número de varreduras até chegar lá muda. Para regras não-monótonas, -o resultado final PODE depender da ordem — cabe a quem escreve a regra -garantir monotonicidade se quiser esse comportamento bem definido. +Repeated sweeps until convergence — for problems with UNBOUNDED SPATIAL +DEPENDENCY (connectivity, flow routing, watershed delineation), where no +fixed-size halo is enough on its own: a cell's value may, in principle, +depend on the whole domain. + +Only the orchestration pattern (small halo + repetition, instead of one +large halo) lives here; the actual rule (e.g. a connectivity propagation +with scipy.ndimage.binary_propagation per block) is supplied by the +caller. Any `rule` that operates on a window with a halo and returns an +updated core works. + +Why this is the right pattern for unbounded dependency +------------------------------------------------------ +A fixed-size halo lets information cross only ONE block boundary per +call. Repeating the sweep N times (N = number of blocks, in the worst +case) lets the propagation "front" advance one block per round — after +N rounds it has necessarily reached every reachable cell in the whole +domain, whatever its size. It is the same principle as a distributed +BFS in bounded rounds. + +Gauss-Seidel, not Jacobi +------------------------ +Each block writes its result IMMEDIATELY back to the same read slot +(through write_block_core_in_place, no ping-pong) — a block processed +later in the SAME sweep already sees the update of a block processed +earlier. This speeds up convergence (fewer sweeps) without changing the +final result, as long as the rule is a monotone operator (e.g. +connectivity only grows, never shrinks) — then the final fixed point +does not depend on block order, only the number of sweeps to reach it +does. For non-monotone rules the final result MAY depend on the order — +it is up to the rule's author to ensure monotonicity if that behaviour +must be well defined. """ from __future__ import annotations @@ -50,38 +47,35 @@ def sweep_until_convergence( max_sweeps: int | None = None, ) -> dict: """ - Repete varreduras de blocos+halo até uma varredura inteira não - mudar nenhum bloco. + Repeat block+halo sweeps until a whole sweep changes no block. Parameters ---------- rule : callable `rule(window: dict[str, np.ndarray]) -> dict[str, np.ndarray]`. - Recebe a janela com halo de cada array do workspace (mesmo - formato de MemmapRasterWorkspace.read_block_with_halo) e - devolve os arrays atualizados, do tamanho do NÚCLEO do bloco - (sem halo) — mesmo contrato das regras usadas em + Receives the halo window of each workspace array (same format as + MemmapRasterWorkspace.read_block_with_halo) and returns the + updated arrays, sized to the block CORE (no halo) — the same + contract as the rules used in HaloChunkedSyncRasterModel/DiskChunkedSyncRasterModel. - boundary_value : escalar ou dict, opcional - Mesmo mecanismo de haloexec.resolve_boundary_value. - max_sweeps : int, opcional - Limite de segurança. Default: número de blocos + 1 (mesmo - limite conservador do protótipo original — cada varredura - completa avança a frente de propagação em pelo menos um - bloco, então esse limite basta para qualquer domínio). + boundary_value : scalar or dict, optional + Same mechanism as haloexec.resolve_boundary_value. + max_sweeps : int, optional + Safety limit. Default: number of blocks + 1 (each complete sweep + advances the propagation front by at least one block, so this + limit is enough for any domain). Returns ------- - dict com "sweeps" (quantas varreduras rodaram até convergir), - "blocks_changed_total" (soma de blocos que mudaram em cada - varredura, ao longo de todas as varreduras). + dict with "sweeps" (how many sweeps ran until convergence), + "blocks_changed_total" (number of blocks changed in each sweep, + summed over all sweeps) and "converged". Raises ------ RuntimeError - Se não convergir dentro de max_sweeps — falha alto em vez de - truncar silenciosamente (mesma disciplina do protótipo - original). + If it does not converge within max_sweeps — fails loudly instead + of truncating silently. """ blocks: list[Block] = workspace.blocks() if max_sweeps is None: @@ -109,7 +103,7 @@ def sweep_until_convergence( return {"sweeps": sweep, "blocks_changed_total": total_changed_blocks, "converged": True} raise RuntimeError( - f"sweep_until_convergence não convergiu em {max_sweeps} varreduras " - f"(limite = número de blocos + 1). Verifique se a regra é monótona " - f"e termina, ou aumente max_sweeps explicitamente." + f"sweep_until_convergence did not converge in {max_sweeps} sweeps " + f"(limit = number of blocks + 1). Check that the rule is monotone " + f"and terminates, or raise max_sweeps explicitly." ) diff --git a/src/haloexec/disk/io/geotiff.py b/src/haloexec/disk/io/geotiff.py index 606f093..2e37d6e 100644 --- a/src/haloexec/disk/io/geotiff.py +++ b/src/haloexec/disk/io/geotiff.py @@ -1,25 +1,21 @@ """ -Carrega um GeoTIFF diretamente para um MemmapRasterWorkspace, bloco a -bloco, via rasterio.windows.Window — nunca materializa uma banda -inteira em RAM. - -Extraído e generalizado do mesmo padrão usado em um protótipo de aluno -(chunked_engine.py, br_mangue_preprocess), que já lia rasters reais -bloco a bloco com `rasterio.open(...).read(banda, window=Window(...))`. -Aqui a extração usa a convenção `band_spec` já estabelecida em -`dissmodel.io.raster.load_geotiff` — lista de (nome, dtype, nodata) — -em vez de nomes de banda hardcoded ("papel", "elevação"), então serve -tanto para o TIFF_BANDS do BR-MANGUE quanto para qualquer outro layout. - -Por que não usar dissmodel.io.raster.load_geotiff diretamente ------------------------------------------------------------------ -Essa função (real, do pacote dissmodel) lê cada banda inteira de uma -vez (`ds.read(i)`) para dentro de um RasterBackend em RAM — correta -para o caminho HaloChunkedSyncRasterModel (que já materializa a grade -global de qualquer forma), mas inadequada para MemmapRasterWorkspace, -cujo propósito é justamente nunca materializar a grade inteira. - -Requer rasterio (extra opcional "geotiff": pip install -e ".[geotiff]"). +Load a GeoTIFF directly into a MemmapRasterWorkspace, block by block, +through rasterio.windows.Window — a whole band is never materialized in +RAM. + +Bands are described with the `band_spec` convention already used by +`dissmodel.io.raster.load_geotiff` — a list of (name, dtype, nodata) — +instead of hard-coded band names, so any band layout works. + +Why not use dissmodel.io.raster.load_geotiff directly +----------------------------------------------------- +That function (from the dissmodel package) reads each whole band at +once (`ds.read(i)`) into an in-RAM RasterBackend — right for the +HaloChunkedSyncRasterModel path (which materializes the global grid +anyway), but unsuitable for MemmapRasterWorkspace, whose whole purpose +is never to materialize the entire grid. + +Requires rasterio (optional "geotiff" extra: pip install -e ".[geotiff]"). """ from __future__ import annotations @@ -45,29 +41,28 @@ def load_geotiff_into_workspace( band_spec: list[tuple[str, str, float]], ) -> None: """ - Popula um MemmapRasterWorkspace bloco a bloco a partir de um único - GeoTIFF. Atalho de conveniência sobre load_geotiffs_into_workspace - (múltiplos arquivos) para o caso comum de um arquivo só. + Fill a MemmapRasterWorkspace block by block from a single GeoTIFF. + Convenience shortcut over load_geotiffs_into_workspace (several + files) for the common one-file case. Parameters ---------- workspace : MemmapRasterWorkspace - Já criado com o mesmo shape do GeoTIFF (workspace.shape deve - bater com (altura, largura) do arquivo) e com os arrays de - `band_spec` já declarados em `arrays=` no `.create()`. + Already created with the GeoTIFF's shape (workspace.shape must + match the file's (height, width)) and with the `band_spec` arrays + declared in `arrays=` of `.create()`. path : str | Path - Caminho do GeoTIFF local. - band_spec : list[(nome, dtype, nodata)] - Mesmo formato de dissmodel.io.raster.load_geotiff — banda 1 - do arquivo mapeia para band_spec[0], banda 2 para band_spec[1], - etc. Arrays cujo nome não foi declarado no workspace são - ignorados (permite carregar só um subconjunto das bandas). - - Nota sobre nodata: este loader NÃO filtra nem substitui valores de - nodata — copia os valores brutos do arquivo. Use o `nodata` de - cada entrada de band_spec para configurar o `boundary_value` - correspondente nas camadas de halo (ver achado documentado no - README sobre boundary_value por-array). + Path of the local GeoTIFF. + band_spec : list[(name, dtype, nodata)] + Same format as dissmodel.io.raster.load_geotiff — file band 1 + maps to band_spec[0], band 2 to band_spec[1], and so on. Arrays + whose name is not declared in the workspace are skipped (so a + subset of the bands can be loaded). + + Note on nodata: this loader does NOT filter or replace nodata + values — it copies the file's raw values. Use each band_spec entry's + `nodata` to set the matching `boundary_value` for the halo layers + (see the README finding on per-array boundary_value). """ load_geotiffs_into_workspace(workspace, [(path, band_spec)]) @@ -77,38 +72,36 @@ def load_geotiffs_into_workspace( sources: list[tuple[str | Path, list[tuple[str, str, float]]]], ) -> None: """ - Popula um MemmapRasterWorkspace bloco a bloco a partir de - MÚLTIPLOS arquivos GeoTIFF, lendo a MESMA janela de bloco de cada - arquivo por vez. - - Generaliza load_geotiff_into_workspace (um único arquivo) para o - padrão de múltiplos rasters usado em um protótipo de aluno - (chunked_engine.py: `dominio_path` + `base_path` lidos juntos, - validando shape/CRS consistentes antes de processar). O motor de - blocos+halo (Block/make_blocks/halo_window) sempre foi agnóstico a - quantos arquivos alimentam a grade — um bloco é só uma posição - (r0:r1, c0:c1); esta função é o que estava faltando para o - carregador acompanhar essa generalidade. + Fill a MemmapRasterWorkspace block by block from SEVERAL GeoTIFF + files, reading the SAME block window from each file at a time. + + Generalizes load_geotiff_into_workspace (a single file) to inputs + split across several rasters that share one grid (e.g. a domain mask + and a base layer read together, with consistent shape/CRS checked + before processing). The block+halo engine (Block/make_blocks/ + halo_window) has always been agnostic to how many files feed the + grid — a block is just a position (r0:r1, c0:c1); this function lets + the loader match that generality. Parameters ---------- workspace : MemmapRasterWorkspace - Já criado com shape batendo com TODOS os arquivos de entrada. - sources : list[(caminho, band_spec)] - Cada arquivo contribui com os arrays declarados em seu próprio - band_spec (mesmo formato de load_geotiff_into_workspace). Um - array pode vir de qualquer um dos arquivos — não precisa haver - sobreposição de nomes entre band_specs de arquivos diferentes. + Already created with a shape matching ALL input files. + sources : list[(path, band_spec)] + Each file contributes the arrays declared in its own band_spec + (same format as load_geotiff_into_workspace). An array may come + from any of the files — band_specs of different files do not + need overlapping names. Raises ------ ValueError - Se os arquivos não tiverem o mesmo shape entre si, ou não - baterem com o shape do workspace, ou tiverem CRS diferentes - entre si (quando CRS está definido em mais de um arquivo). + If the files do not share the same shape, do not match the + workspace shape, or have different CRSs (when a CRS is defined in + more than one file). """ if not HAS_RASTERIO: - raise ImportError("rasterio é necessário — pip install -e '.[geotiff]'") + raise ImportError("rasterio is required — pip install -e '.[geotiff]'") declared = set(workspace.metadata["arrays"]) datasets = [rasterio.open(str(path)) for path, _ in sources] @@ -120,19 +113,19 @@ def load_geotiffs_into_workspace( shape = (ds.height, ds.width) if shape != ref_shape: raise ValueError( - f"Shape inconsistente entre arquivos: {path} tem {shape}, " - f"esperado {ref_shape} (do primeiro arquivo da lista)." + f"Inconsistent shape between files: {path} has {shape}, " + f"expected {ref_shape} (from the first file in the list)." ) if ds.crs is not None and ref_crs is not None and ds.crs != ref_crs: raise ValueError( - f"CRS inconsistente entre arquivos: {path} tem {ds.crs}, " - f"esperado {ref_crs} (do primeiro arquivo da lista)." + f"Inconsistent CRS between files: {path} has {ds.crs}, " + f"expected {ref_crs} (from the first file in the list)." ) if ref_shape != tuple(workspace.shape): raise ValueError( - f"Shape dos GeoTIFFs {ref_shape} não bate com o shape " - f"do workspace {tuple(workspace.shape)}." + f"Shape of the GeoTIFFs {ref_shape} does not match the " + f"workspace shape {tuple(workspace.shape)}." ) for block in workspace.blocks(): @@ -159,38 +152,41 @@ def save_workspace_to_geotiff( compress: str = "lzw", ) -> None: """ - Grava arrays do MemmapRasterWorkspace em um GeoTIFF em janelas - (bloco a bloco), sem nunca materializar a grade inteira em RAM. + Write MemmapRasterWorkspace arrays to a GeoTIFF in windows (block + by block), never materializing the whole grid in RAM. Parameters ---------- workspace : MemmapRasterWorkspace - Workspace de onde ler os dados do slot de leitura atual. + Workspace whose current read slot is written. path : str | Path - Caminho do arquivo GeoTIFF de saída. - bands : list[str] ou list[(nome, dtype, nodata)] - Nomes dos arrays a gravar como bandas (1, 2, ...). + Path of the output GeoTIFF. + bands : list[str] or list[(name, dtype, nodata)] + Names of the arrays to write as bands (1, 2, ...). transform : Affine, optional - Matriz de geotransformação do rasterio. Se None, cria uma padrão. - crs : CRS ou str, default="EPSG:31984" - Sistema de referência de coordenadas. + rasterio geotransform. If None, a default one is used (with a + warning). + crs : CRS or str, default="EPSG:31984" + Coordinate reference system. compress : str, default="lzw" - Compressão do GeoTIFF. + GeoTIFF compression. + Note ---- - GeoTIFF/GDAL exige um único dtype e um único nodata para TODAS as - bandas de um arquivo (limitação do formato, não deste código — - testado empiricamente: GDAL colapsa nodata por banda para um valor - só, mesmo passando valores distintos). Se `bands` misturar dtypes - ou nodata diferentes, esta função converte tudo para o dtype comum - (menor tipo que comporta todos, via `np.result_type`) e usa o - nodata da PRIMEIRA banda para o arquivo inteiro — e emite um aviso - (`warnings.warn`) sempre que isso descartar informação, em vez de - fazer silenciosamente. Se cada array precisa manter seu próprio - dtype/nodata, grave um GeoTIFF por array em vez de multibanda. + GeoTIFF/GDAL requires a single dtype and a single nodata for ALL + bands of a file (a limitation of the format, not of this code — + tested empirically: GDAL collapses per-band nodata to one value, + even when distinct values are passed). If `bands` mixes dtypes or + nodata values, this function converts everything to the common dtype + (the smallest type that holds them all, via `np.result_type`) and + uses the FIRST band's nodata for the whole file — and emits a + warning (`warnings.warn`) whenever that discards information, rather + than doing it silently. If each array must keep its own + dtype/nodata, write one GeoTIFF per array instead of a multi-band + file. """ if not HAS_RASTERIO: - raise ImportError("rasterio é necessário — pip install -e '.[geotiff]'") + raise ImportError("rasterio is required — pip install -e '.[geotiff]'") import warnings @@ -212,11 +208,11 @@ def save_workspace_to_geotiff( if transform is None: transform = from_origin(500_000.0, 9_700_000.0, 30.0, 30.0) warnings.warn( - "save_workspace_to_geotiff: nenhum `transform` foi passado — " - "usando origem padrão arbitrária (500000, 9700000, 30x30m). " - "O GeoTIFF resultante NÃO estará georreferenciado corretamente " - "a menos que essa origem coincida com o domínio real. Passe " - "`transform=` explicitamente para dados reais.", + "save_workspace_to_geotiff: no `transform` was given — " + "using an arbitrary default origin (500000, 9700000, 30x30 m). " + "The resulting GeoTIFF will NOT be correctly georeferenced " + "unless that origin matches the real domain. Pass " + "`transform=` explicitly for real data.", stacklevel=2, ) @@ -224,12 +220,12 @@ def save_workspace_to_geotiff( common_dtype = np.result_type(*(np.dtype(d) for d in distinct_dtypes)) if len(distinct_dtypes) > 1: warnings.warn( - f"save_workspace_to_geotiff: bandas com dtypes distintos " - f"{sorted(distinct_dtypes)} — GeoTIFF exige um único dtype por " - f"arquivo, convertendo tudo para {common_dtype} (pode aumentar " - f"o tamanho do arquivo e/ou alterar a semântica de arrays " - f"categóricos). Para preservar dtypes, grave um GeoTIFF " - f"separado por array.", + f"save_workspace_to_geotiff: bands with different dtypes " + f"{sorted(distinct_dtypes)} — GeoTIFF requires a single dtype per " + f"file, converting everything to {common_dtype} (may increase " + f"the file size and/or change the meaning of categorical " + f"arrays). To keep the dtypes, write a separate GeoTIFF " + f"per array.", stacklevel=2, ) @@ -237,12 +233,12 @@ def save_workspace_to_geotiff( first_nodata = parsed_bands[0][2] if len(distinct_nodatas) > 1: warnings.warn( - f"save_workspace_to_geotiff: bandas com nodata distintos " - f"{sorted(distinct_nodatas)} — GeoTIFF só suporta um nodata " - f"por arquivo (limitação do GDAL, testado empiricamente: " - f"valores por banda são silenciosamente colapsados). Usando " - f"nodata da primeira banda ({first_nodata!r}) para o arquivo " - f"inteiro; as demais bandas ficam sem nodata correto.", + f"save_workspace_to_geotiff: bands with different nodata values " + f"{sorted(distinct_nodatas)} — GeoTIFF supports only one nodata " + f"per file (a GDAL limitation, tested empirically: per-band " + f"values are silently collapsed). Using the first band's " + f"nodata ({first_nodata!r}) for the whole file; the other " + f"bands will not have the correct nodata.", stacklevel=2, ) diff --git a/src/haloexec/disk/io/zarr.py b/src/haloexec/disk/io/zarr.py index 982d99f..66bcece 100644 --- a/src/haloexec/disk/io/zarr.py +++ b/src/haloexec/disk/io/zarr.py @@ -1,28 +1,28 @@ """ -Carrega arrays de um Zarr store direto para MemmapRasterWorkspace, -bloco a bloco — segunda opção de entrada de dados, ao lado de -geotiff_io.py (TIFF/VRT). - -Motivação: o disscube (DisSModel/disscube) armazena variáveis -derivadas nativamente em Zarr (`data/derived/{grid_id}/{tile_id}/ -{spec_hash}/{variable_name}.zarr`), já alinhadas à grade mestra pelo -seu GridAligner (com resampling por-operador e alinhamento fino para -categóricas — ver README). Este módulo permite consumir esse dado -direto, sem precisar materializar para GeoTIFF primeiro. - -Por que um módulo separado de geotiff_io.py ---------------------------------------------- -Zarr já é nativamente chunked — não precisa de rasterio.windows.Window -nem de VRT para leitura parcial; um zarr.Array suporta slicing direto -(`arr[r0:r1, c0:c1]`) que só lê os chunks necessários do store. A API -pública espelha geotiff_io.py (mesma forma de declarar quais arrays do -workspace vêm de qual variável de origem), para que os dois caminhos -de entrada (GeoTIFF/VRT e Zarr) sejam intercambiáveis do ponto de -vista de quem usa MemmapRasterWorkspace — troca-se o loader, o resto -do pipeline (halo, disco, modelos) não muda. - -Requer o extra opcional "zarr" (zarr>=2.16). xarray/rioxarray NÃO são -necessários aqui — lê-se o zarr.Array bruto por nome de variável. +Load arrays from a Zarr store directly into a MemmapRasterWorkspace, +block by block — the second data-input path, next to geotiff.py +(GeoTIFF/VRT). + +Motivation: disscube (DisSModel/disscube) stores derived variables +natively in Zarr (`data/derived/{grid_id}/{tile_id}/{spec_hash}/ +{variable_name}.zarr`), already aligned to the master grid by its +GridAligner (per-operator resampling and fine alignment for categorical +data — see its README). This module consumes that data directly, +without materializing it to GeoTIFF first. + +Why a module separate from geotiff.py +------------------------------------- +Zarr is natively chunked — it needs neither rasterio.windows.Window +nor a VRT for partial reads; a zarr.Array supports direct slicing +(`arr[r0:r1, c0:c1]`) that reads only the chunks it needs from the +store. The public API mirrors geotiff.py (same way of declaring which +workspace arrays come from which source variable), so the two input +paths (GeoTIFF/VRT and Zarr) are interchangeable for code that uses +MemmapRasterWorkspace — swap the loader, and the rest of the pipeline +(halo, disk, models) stays the same. + +Requires the optional "zarr" extra (zarr>=3). xarray/rioxarray are NOT +needed here — the raw zarr.Array is read by variable name. """ from __future__ import annotations @@ -41,26 +41,25 @@ def _resolve_axis_order(arr, expected_names: tuple[str, ...]) -> tuple[int, ...] | None: - """Usa arr.metadata.dimension_names (campo nativo do Zarr v3, o - mesmo que xarray grava ao salvar um DataArray via to_zarr) para - determinar a ordem real dos eixos no array em disco, e devolve os - índices de transposição necessários para chegar em expected_names. - - Por que isso é necessário: xarray/disscube NÃO garantem que a - ordem de eixos gravada em disco seja (y, x) — o próprio - CubeClient.load() do disscube faz `.transpose("y", "x")` - defensivamente antes de usar qualquer array, precisamente porque - a ordem pode vir diferente. Um array QUADRADO com eixos trocados - tem o MESMO shape nos dois casos — a checagem de shape sozinha - não detecta a troca; é silenciosa, não trava. - - No Zarr v2 (zarr-python 2.x, o único disponível em Python 3.10) - não existe dimension_names: o xarray grava os nomes no atributo - `_ARRAY_DIMENSIONS`, que é lido como fallback. - - Retorna None se nenhum dos dois estiver disponível (Zarr sem - metadado de dimensão — não há como verificar, assume-se a ordem - como está, mesmo comportamento de antes desta correção). + """Use arr.metadata.dimension_names (the native Zarr v3 field, the + same one xarray writes when saving a DataArray with to_zarr) to find + the actual axis order of the array on disk, and return the transpose + indices that bring it to expected_names. + + Why this is needed: xarray/disscube do NOT guarantee that the axis + order written to disk is (y, x) — disscube's own CubeClient.load() + calls `.transpose("y", "x")` defensively before using any array, + precisely because the order may differ. A SQUARE array with swapped + axes has the SAME shape either way — a shape check alone does not + catch the swap; it is silent, nothing fails. + + Zarr v2-format arrays have no dimension_names: xarray stores the + names in the `_ARRAY_DIMENSIONS` attribute, which is read as a + fallback. + + Returns None if neither is available (a Zarr store without dimension + metadata — there is nothing to check against, so the order is taken + as it is). """ dims = getattr(getattr(arr, "metadata", None), "dimension_names", None) if not dims: @@ -70,23 +69,23 @@ def _resolve_axis_order(arr, expected_names: tuple[str, ...]) -> tuple[int, ...] return None dims = tuple(dims) if set(dims) != set(expected_names): - return None # nomes de dimensão inesperados -- não arrisca reordenar + return None # unexpected dimension names -- do not risk reordering return tuple(dims.index(name) for name in expected_names) def _open_variable(store: str, variable_name: str | None): - """Abre um zarr store, que pode ser um grupo (com várias variáveis, - acessadas por nome) ou um array único (variável direta).""" + """Open a Zarr store, which may be a group (several variables, + accessed by name) or a single array (the variable itself).""" opened = zarr.open(store, mode="r") if hasattr(opened, "arrays") or hasattr(opened, "array_keys"): - # É um grupo (zarr.Group) — precisa do nome da variável dentro dele. + # A group (zarr.Group) — needs the name of the variable inside it. if variable_name is None: raise ValueError( - f"'{store}' é um grupo Zarr com múltiplas variáveis — " - f"informe o nome da variável em variable_map." + f"'{store}' is a Zarr group with several variables — " + f"give the variable name in variable_map." ) return opened[variable_name] - # É um array único — variable_name é ignorado (ou usado só como rótulo). + # A single array — variable_name is ignored (or used only as a label). return opened @@ -97,33 +96,33 @@ def load_zarr_into_workspace( time_index: int | None = None, ) -> None: """ - Popula um MemmapRasterWorkspace bloco a bloco a partir de um Zarr - store, lendo apenas a fatia de cada bloco por vez. + Fill a MemmapRasterWorkspace block by block from a Zarr store, + reading only one block's slice at a time. Parameters ---------- store : str | Path - Caminho do Zarr store — pode ser um grupo (múltiplas variáveis, - acessadas por nome) ou um array único. - variable_map : dict[nome_array_workspace, nome_variavel_zarr], optional - Mapeia nomes de array do workspace para nomes de variável - dentro do grupo zarr. Se None, assume que os nomes já batem - (mesmo nome no workspace e no zarr), e que `store` é um único - array (não um grupo) se nenhuma variável for nomeada. + Path to the Zarr store — either a group (several variables, + accessed by name) or a single array. + variable_map : dict[workspace_array_name, zarr_variable_name], optional + Maps workspace array names to variable names inside the Zarr + group. If None, the names are assumed to match (same name in the + workspace and in the store), and `store` is taken as a single + array (not a group) when no variable is named. time_index : int, optional - Se a variável tiver uma dimensão temporal inicial (shape - (time, y, x), padrão do "Temporal Backend" do disscube para - produtos derivados com janela de validade), qual índice de - tempo carregar. Obrigatório se a variável for 3D. + If the variable has a leading time dimension (shape + (time, y, x), the layout of disscube's "Temporal Backend" for + derived products with a validity window), which time index to + load. Required when the variable is 3-D. Raises ------ ValueError - Se o shape (y, x) da variável não bater com workspace.shape, - ou se uma variável 3D for informada sem time_index. + If the variable's (y, x) shape does not match workspace.shape, + or if a 3-D variable is given without time_index. """ if not HAS_ZARR: - raise ImportError("zarr é necessário — pip install -e '.[zarr]'") + raise ImportError("zarr is required — pip install -e '.[zarr]'") store = str(store) declared = set(workspace.metadata["arrays"]) @@ -138,23 +137,22 @@ def load_zarr_into_workspace( if arr.ndim == 3: if time_index is None: raise ValueError( - f"Variável '{zarr_var_name}' tem 3 dimensões (provável " - f"dimensão temporal) — informe time_index." + f"Variable '{zarr_var_name}' has 3 dimensions (probably " + f"a time dimension) — give time_index." ) expected_dims = ("time", "y", "x") elif arr.ndim == 2: expected_dims = ("y", "x") else: raise ValueError( - f"Variável '{zarr_var_name}' tem {arr.ndim} dimensões; esperado 2 ou 3." + f"Variable '{zarr_var_name}' has {arr.ndim} dimensions; expected 2 or 3." ) - # Normaliza a ordem de eixos para (y, x) ou (time, y, x) usando - # o metadado nativo de dimensão do Zarr v3, quando disponível. - # Necessário porque a ordem de eixos gravada NÃO é garantida - # (ver docstring de _resolve_axis_order) — um array quadrado - # com eixos trocados tem o mesmo shape nos dois casos, então a - # checagem de shape sozinha não pega a inversão. + # Normalize the axis order to (y, x) or (time, y, x) using the + # store's dimension metadata, when available. Needed because the + # written axis order is NOT guaranteed (see the docstring of + # _resolve_axis_order) — a square array with swapped axes has the + # same shape either way, so a shape check alone misses the swap. axis_order = _resolve_axis_order(arr, expected_dims) shape2d = arr.shape[1:] if arr.ndim == 3 else arr.shape @@ -164,8 +162,8 @@ def load_zarr_into_workspace( if shape2d != tuple(workspace.shape): raise ValueError( - f"Shape de '{zarr_var_name}' {shape2d} não bate com o " - f"shape do workspace {tuple(workspace.shape)}." + f"Shape of '{zarr_var_name}' {shape2d} does not match the " + f"workspace shape {tuple(workspace.shape)}." ) for block in workspace.blocks(): @@ -188,9 +186,9 @@ def load_zarr_into_workspace( raw = np.asarray(arr[tuple(disk_index)]) if axis_order is not None and raw.ndim == 2: - # raw ainda está na ordem relativa em disco (menos o - # eixo de tempo, já reduzido pela indexação inteira - # acima) -- transpõe para (y, x) canônico. + # raw is still in the on-disk relative order (minus the + # time axis, already dropped by the integer indexing + # above) -- transpose to canonical (y, x). def _shift(pos, time_disk_axis=time_disk_axis): return pos - 1 if (time_disk_axis is not None and time_disk_axis < pos) else pos data = np.transpose(raw, (_shift(y_disk_axis), _shift(x_disk_axis))) @@ -210,120 +208,120 @@ def load_zarr_tiles_into_workspace( skip_empty_blocks: bool = False, ) -> None: """ - Popula UM array do workspace a partir de N stores Zarr posicionados - lado a lado — o caso multi-tile, que `load_zarr_into_workspace` não - cobre (ela recebe um store e exige que ele tenha o shape do workspace - inteiro). + Fill ONE workspace array from N Zarr stores placed side by side — + the multi-tile case, which `load_zarr_into_workspace` does not cover + (it takes one store and requires it to have the whole workspace + shape). - Está para o Zarr como o VRT do geomosaic está para o GeoTIFF: junta - pedaços numa grade contínua. A diferença é que não existe formato de - mosaico para Zarr, então a costura acontece aqui, na leitura. + It is to Zarr what geomosaic's VRT is to GeoTIFF: it joins pieces + into one continuous grid. The difference is that Zarr has no mosaic + format, so the stitching happens here, at read time. Parameters ---------- tiles : list[dict] - Um dicionário por pedaço, com as chaves: - - - ``url`` — caminho do store Zarr - - ``variable`` — nome da variável dentro do store - - ``row_off`` — linha, em pixel, onde o pedaço começa no workspace - - ``col_off`` — coluna, em pixel - - ``height`` — altura do pedaço, em pixel - - ``width`` — largura do pedaço, em pixel - - É exatamente o formato que ``CubeClient.tile_layout()`` do - disscube devolve, mas nada aqui depende do disscube: qualquer - origem que saiba dizer caminho e posição serve. Chaves extras são - ignoradas. + One dictionary per piece, with the keys: + + - ``url`` — path of the Zarr store + - ``variable`` — name of the variable inside the store + - ``row_off`` — row, in pixels, where the piece starts in the workspace + - ``col_off`` — column, in pixels + - ``height`` — piece height, in pixels + - ``width`` — piece width, in pixels + + This is exactly the format returned by disscube's + ``CubeClient.tile_layout()``, but nothing here depends on + disscube: any source that can give a path and a position works. + Extra keys are ignored. array : str, optional - Nome do array NO WORKSPACE a preencher. Se None, usa o - ``variable`` do primeiro tile — útil quando os nomes coincidem. + Name of the WORKSPACE array to fill. If None, the ``variable`` of + the first tile is used — handy when the names match. fill : float, optional - Valor para as células que nenhum tile cobre (buracos da malha, - cantos fora da área de estudo). Se None, usa NaN para arrays de - ponto flutuante e 0 para inteiros. + Value for cells no tile covers (holes in the tiling, corners + outside the study area). If None, NaN for floating-point arrays + and 0 for integer ones. skip_empty_blocks : bool - Se True, blocos que nenhum tile toca não são escritos. Como os - `.dat` do workspace nascem esparsos, isso deixa esses blocos sem - ocupar disco — mas eles passam a LER COMO ZERO, não como `fill`. - Só use quando 0 não for um valor válido do domínio (ver a nota - "Esparsidade e custo de disco" no README). Default False, que - escreve `fill` e mantém a distinção ao custo do disco. + If True, blocks that no tile touches are not written. Since the + workspace `.dat` files start sparse, those blocks then take no + disk space — but they READ AS ZERO, not as `fill`. Use only when + 0 is not a valid value in the domain (see "Sparsity and disk + cost" in the README). Default False, which writes `fill` and + keeps the distinction at the cost of disk space. Raises ------ ValueError - Se `tiles` estiver vazio, se o array não for declarado no - workspace, se algum tile faltar chave obrigatória, ou se um tile - cair fora dos limites do workspace — todos casos em que seguir - adiante produziria um mosaico silenciosamente errado. + If `tiles` is empty, if the array is not declared in the + workspace, if a tile lacks a required key, or if a tile falls + outside the workspace bounds — every case in which carrying on + would produce a silently wrong mosaic. ImportError - Se o extra "zarr" não estiver instalado. + If the "zarr" extra is not installed. """ if not HAS_ZARR: - raise ImportError("zarr é necessário — pip install -e '.[zarr]'") + raise ImportError("zarr is required — pip install -e '.[zarr]'") if not tiles: - raise ValueError("A lista de tiles está vazia — nada a carregar.") + raise ValueError("The list of tiles is empty — nothing to load.") - obrigatorias = {"url", "variable", "row_off", "col_off", "height", "width"} + required = {"url", "variable", "row_off", "col_off", "height", "width"} for i, t in enumerate(tiles): - faltando = obrigatorias - set(t) - if faltando: + missing = required - set(t) + if missing: raise ValueError( - f"tiles[{i}] não tem as chaves {sorted(faltando)}; " - f"cada tile precisa de {sorted(obrigatorias)}." + f"tiles[{i}] is missing the keys {sorted(missing)}; " + f"every tile needs {sorted(required)}." ) - # Dois pedaços na mesma posição não é ambiguidade a resolver por ordem: - # um sobrescreveria o outro em silêncio. Acontece de verdade quando o - # layout mistura fatias temporais da mesma variável — cada ano repete as - # mesmas posições. - ocupadas: dict[tuple[int, int], str] = {} + # Two pieces at the same position are not an ambiguity to resolve by + # order: one would silently overwrite the other. It really happens when + # the layout mixes time slices of the same variable — each year repeats + # the same positions. + occupied: dict[tuple[int, int], str] = {} for t in tiles: - chave = (t["row_off"], t["col_off"]) - if chave in ocupadas: + key = (t["row_off"], t["col_off"]) + if key in occupied: raise ValueError( - f"Dois tiles ocupam a posição ({chave[0]}, {chave[1]}): " - f"{ocupadas[chave]!r} e {t.get('tile_id')!r}. Um sobrescreveria " - f"o outro. Se o layout mistura fatias temporais, escolha uma " - f"antes de carregar." + f"Two tiles occupy position ({key[0]}, {key[1]}): " + f"{occupied[key]!r} and {t.get('tile_id')!r}. One would overwrite " + f"the other. If the layout mixes time slices, pick one " + f"before loading." ) - ocupadas[chave] = t.get("tile_id") + occupied[key] = t.get("tile_id") array_name = array or tiles[0]["variable"] - declarados = set(workspace.metadata["arrays"]) - if array_name not in declarados: + declared = set(workspace.metadata["arrays"]) + if array_name not in declared: raise ValueError( - f"O workspace não declara o array {array_name!r} " - f"(declarados: {sorted(declarados)})." + f"The workspace does not declare the array {array_name!r} " + f"(declared: {sorted(declared)})." ) - altura, largura = workspace.shape + height, width = workspace.shape for t in tiles: if (t["row_off"] < 0 or t["col_off"] < 0 - or t["row_off"] + t["height"] > altura - or t["col_off"] + t["width"] > largura): + or t["row_off"] + t["height"] > height + or t["col_off"] + t["width"] > width): raise ValueError( - f"tile {t.get('tile_id')!r} em " + f"tile {t.get('tile_id')!r} at " f"({t['row_off']},{t['col_off']}) {t['height']}x{t['width']} " - f"não cabe no workspace {altura}x{largura}." + f"does not fit in the workspace {height}x{width}." ) dtype = np.dtype(workspace.metadata["arrays"][array_name]) if fill is None: fill = np.nan if np.issubdtype(dtype, np.floating) else 0 - abertos = {} + opened = {} try: for t in tiles: - chave = (t["url"], t["variable"]) - if chave not in abertos: - abertos[chave] = _open_variable(t["url"], t["variable"]) - - # Percorre os BLOCOS do workspace, não os tiles: um bloco pode cair - # sobre dois tiles vizinhos, ou sobre um buraco da malha. Montá-lo a - # partir de tudo que o cobre é o que faz a costura ficar correta — - # escrever tile a tile deixaria as bordas dependendo da ordem. + key = (t["url"], t["variable"]) + if key not in opened: + opened[key] = _open_variable(t["url"], t["variable"]) + + # Iterate over the workspace BLOCKS, not the tiles: a block may span + # two neighbouring tiles, or a hole in the tiling. Assembling it from + # everything that covers it is what makes the stitching correct — + # writing tile by tile would make the edges depend on the order. for block in workspace.blocks(): buf = None for t in tiles: @@ -339,17 +337,17 @@ def load_zarr_tiles_into_workspace( (block.r1 - block.r0, block.c1 - block.c0), fill, dtype=dtype ) - arr = abertos[(t["url"], t["variable"])] - ordem = _resolve_axis_order(arr, ("y", "x")) - trecho = np.asarray(arr[ + arr = opened[(t["url"], t["variable"])] + order = _resolve_axis_order(arr, ("y", "x")) + piece = np.asarray(arr[ r0 - t["row_off"]:r1 - t["row_off"], c0 - t["col_off"]:c1 - t["col_off"], - ]) if ordem in (None, (0, 1)) else np.asarray(arr[ + ]) if order in (None, (0, 1)) else np.asarray(arr[ c0 - t["col_off"]:c1 - t["col_off"], r0 - t["row_off"]:r1 - t["row_off"], ]).T - buf[r0 - block.r0:r1 - block.r0, c0 - block.c0:c1 - block.c0] = trecho + buf[r0 - block.r0:r1 - block.r0, c0 - block.c0:c1 - block.c0] = piece if buf is None: if skip_empty_blocks: @@ -360,6 +358,6 @@ def load_zarr_tiles_into_workspace( workspace.write_block_to_read_slot(block, array_name, buf) finally: - abertos.clear() + opened.clear() workspace.flush() diff --git a/src/haloexec/disk/sync_model.py b/src/haloexec/disk/sync_model.py index 10bbfdd..1e4cc78 100644 --- a/src/haloexec/disk/sync_model.py +++ b/src/haloexec/disk/sync_model.py @@ -1,34 +1,34 @@ """ -Integração disco + halo com SyncRasterModel (FloodModel, MangroveModel, -e qualquer outro modelo dissmodel do mesmo padrão), sem modificar o -pacote dissmodel instalado. - -Por que não basta reaproveitar HaloChunkedSyncRasterModel (in-memory) +Disk + halo integration for SyncRasterModel (FloodModel, MangroveModel, +and any other dissmodel model following the same pattern), without +modifying the installed dissmodel package. + +Why reusing HaloChunkedSyncRasterModel (in-memory) is not enough +---------------------------------------------------------------- +SyncRasterModel.pre_execute()/post_execute() call synchronize(), which +does `self.backend.get(name).copy()` on the WHOLE array — if +`self.backend` were a wrapper over an np.memmap, that `.copy()` would +materialize the whole grid in RAM just to take the "_past" +snapshot, defeating the purpose of using disk. + +This module REPLICATES SyncRasterModel's "_past" synchronization logic +locally (it neither imports nor modifies it in the installed package), +adapted to copy block by block between memmaps through +MemmapRasterWorkspace. When the migration to the dissmodel core +happens, this is the part to reconcile with +dissmodel.geo.raster.sync_model.SyncRasterModel.synchronize() — for +now the two live side by side, with no coupling. + +Double-buffer semantics for "_past" +----------------------------------- +Unlike the "current" arrays (e.g. "uso", "alt"), which only become +ready in the OTHER slot after a complete step (classic ping-pong), the +"_past" arrays must be available in the SAME slot that execute() will +READ in that step — which is why they use write_block_to_read_slot(), +not write_block_core(). See the MemmapRasterWorkspace docstring. + +Usage (array and parameter names are those of the BR-MANGUE FloodModel) ----------------------------------------------------------------------- -SyncRasterModel.pre_execute()/post_execute() chamam synchronize(), que -faz `self.backend.get(name).copy()` sobre o array INTEIRO — se -`self.backend` fosse um wrapper em cima de um np.memmap, esse `.copy()` -materializaria a grade inteira em RAM só para gerar o snapshot -"_past", anulando o propósito de usar disco. - -Este módulo REPLICA localmente a lógica de sincronização "_past" do -SyncRasterModel (não a importa nem a modifica no pacote instalado), -adaptada para copiar bloco a bloco entre memmaps via -MemmapRasterWorkspace. Quando a migração para dissmodel core acontecer, -este é o trecho que deve ser reconciliado com -dissmodel.geo.raster.sync_model.SyncRasterModel.synchronize() — por -ora, os dois vivem em paralelo, sem acoplamento. - -Semântica de double-buffer para "_past" ------------------------------------------ -Diferente dos arrays "correntes" (ex. "uso", "alt"), que só ficam -prontos no OUTRO slot após um passo completo (ping-pong clássico), os -arrays "_past" precisam estar disponíveis no MESMO slot que execute() -vai LER naquele passo — por isso usam write_block_to_read_slot(), não -write_block_core(). Ver docstring de MemmapRasterWorkspace. - -Uso ---- class FloodModelDiskHalo(DiskChunkedSyncRasterModel, FloodModel): pass @@ -58,9 +58,9 @@ def workspace_arrays_for_sync_model( base: dict[str, np.dtype], land_use_types: list[str], ) -> dict[str, np.dtype]: - """Deriva o dict de arrays a declarar em MemmapRasterWorkspace.create(), - adicionando automaticamente "_past" para cada nome em - land_use_types (mesmo dtype do array base).""" + """Build the dict of arrays to declare in MemmapRasterWorkspace.create(), + adding "_past" automatically for each name in land_use_types + (same dtype as the base array).""" arrays = dict(base) for name in land_use_types: arrays[f"{name}_past"] = np.dtype(base[name]) @@ -69,12 +69,11 @@ def workspace_arrays_for_sync_model( class DiskChunkedSyncRasterModel: """ - Mixin que processa um SyncRasterModel (ex.: FloodModel) em blocos - lidos de um MemmapRasterWorkspace, incluindo a sincronização - "_past" feita bloco a bloco — sem nunca materializar a grade - inteira em RAM. + Mixin that runs a SyncRasterModel (e.g. FloodModel) in blocks read + from a MemmapRasterWorkspace, including the "_past" synchronization + done block by block — never materializing the whole grid in RAM. - Ordem de herança (MRO): este mixin deve vir primeiro, ex.: + Inheritance order (MRO): this mixin must come first, e.g. `class FloodModelDiskHalo(DiskChunkedSyncRasterModel, FloodModel)`. """ @@ -85,17 +84,17 @@ def setup(self, workspace: MemmapRasterWorkspace, halo: int | None = None, self.boundary_value = boundary_value self._synced_before_first_execute = False - # Placeholder leve: RasterBackend(shape=...) não aloca arrays, - # só existe para satisfazer o contrato de RasterModel.setup() - # (self.backend = backend; self.shape = backend.shape). O - # backend real por bloco é criado dentro de execute(). + # Lightweight placeholder: RasterBackend(shape=...) allocates no + # arrays; it only satisfies RasterModel.setup()'s contract + # (self.backend = backend; self.shape = backend.shape). The real + # per-block backend is created inside execute(). placeholder = RasterBackend(shape=workspace.shape) - super().setup(backend=placeholder, **kwargs) # delega para a subclasse real + super().setup(backend=placeholder, **kwargs) # delegate to the real subclass def _synchronize_via_workspace(self) -> None: - """Equivalente bloco-a-bloco de SyncRasterModel.synchronize(): - copia "" -> "_past", dentro do MESMO slot de - leitura atual (ver docstring do módulo).""" + """Block-by-block equivalent of SyncRasterModel.synchronize(): + copies "" -> "_past" within the SAME current read + slot (see the module docstring).""" for name in getattr(self, "land_use_types", []): for block in self.workspace.blocks(): values = self.workspace.read_block_core(block, name) @@ -124,22 +123,22 @@ def execute(self) -> None: self.backend = block_backend self.shape = block_backend.shape try: - super().execute() # lógica real (ex.: FloodModel.execute) + super().execute() # the real logic (e.g. FloodModel.execute) finally: self.backend = real_backend self.shape = real_shape updates = {} for name, arr in block_backend.arrays.items(): - # IMPORTANTE: não excluir "_past" aqui. Se este modelo não - # gerencia um determinado "_past" (ex.: FloodModel não - # gerencia "solo_past", só MangroveModel gerencia), ele - # ainda precisa ser levado adiante sem alteração através - # do swap — senão fica órfão no slot novo (nunca escrito, - # permanece com o valor zerado/obsoleto da alocação - # inicial do memmap). O "_past" que ESTE modelo gerencia - # será corretamente sobrescrito por _synchronize_via_workspace - # em post_execute(), já no slot pós-swap. + # IMPORTANT: do not skip "_past" here. If this model does + # not manage a given "_past" (e.g. FloodModel does not + # manage "solo_past", only MangroveModel does), it still + # has to be carried over unchanged across the swap — + # otherwise it is orphaned in the new slot (never written, + # left with the zeroed/stale value of the memmap's initial + # allocation). The "_past" that THIS model manages is + # correctly overwritten by _synchronize_via_workspace in + # post_execute(), already in the post-swap slot. core = arr[h:-h, h:-h] if h > 0 else arr updates[name] = core ws.write_block_core(block, updates) diff --git a/src/haloexec/disk/workspace.py b/src/haloexec/disk/workspace.py index 9cbebd9..645a499 100644 --- a/src/haloexec/disk/workspace.py +++ b/src/haloexec/disk/workspace.py @@ -1,35 +1,32 @@ """ -Workspace de arrays em disco (np.memmap) com decomposição em blocos+halo, -double-buffering e checkpoint — para grades grandes demais para caber -inteiras em RAM. - -Extraído e generalizado de um protótipo de aluno (chunked_engine.py, -br_mangue_preprocess) que tinha a ideia certa mas amarrada a nomes de -estado específicos do domínio BR-MANGUE. Esta versão é genérica: não -sabe nada sobre "papel", "uso", "alt" ou qualquer domínio — apenas -armazena arrays 2D nomeados com dtype declarado. - -Por que existe separado do resto do haloexec ----------------------------------------------- -`HaloChunkedRasterCellularAutomaton` e `HaloChunkedSyncRasterModel` -(dissmodel_ca.py, sync_model.py) fazem halo via `np.pad` da grade -INTEIRA em memória — funciona bem até a grade caber em RAM. Quando não -cabe (o caso real de domínios costa Pará-Maranhão em escala), é preciso -nunca materializar a grade inteira: ler só a janela do bloco+halo do -disco, recortada nas bordas quando o bloco toca a borda do domínio. - -Fundamentação teórica: mesma do resto do haloexec (Kjolstad & Snir -2010; Xia et al. 2025) — aqui aplicada com a variante em que o halo é -lido diretamente do arquivo em disco, sem padding em memória da grade -global. +On-disk array workspace (np.memmap) with block+halo decomposition, +double-buffering and checkpointing — for grids too large to fit in RAM +as a whole. + +Generic: it knows nothing about any model's state variables or domain — +it only stores named 2-D arrays with a declared dtype. + +Why it is separate from the rest of haloexec +-------------------------------------------- +`HaloChunkedRasterCellularAutomaton` and `HaloChunkedSyncRasterModel` +(ram/) build the halo with `np.pad` over the WHOLE grid in memory — +fine while the grid fits in RAM. When it does not (e.g. a stretch of +coastline at fine resolution), the whole grid must never be +materialized: only the block+halo window is read from disk, clipped at +the edges when the block touches the domain boundary. + +Theoretical basis: the same as the rest of haloexec (Kjolstad & Snir +2010; Xia et al. 2025), here in the variant where the halo is read +directly from the file on disk, with no in-memory padding of the global +grid. Double-buffering ------------------ -Cada array tem dois slots físicos ("a" e "b"). Um passo de tempo lê do -slot corrente e escreve no outro slot — evita que uma célula leia o -valor já atualizado de um vizinho no mesmo passo (hazard clássico de -autômatos celulares síncronos). Ao final do passo, os slots trocam de -papel (swap lógico, sem copiar dados). +---------------- +Each array has two physical slots ("a" and "b"). A time step reads from +the current slot and writes to the other one — so a cell never reads a +neighbour's value already updated in the same step (the classic hazard +of synchronous cellular automata). At the end of the step the slots +swap roles (a logical swap, no data copied). """ from __future__ import annotations @@ -53,20 +50,20 @@ def _write_json_atomic(path: Path, payload: dict[str, Any]) -> None: @dataclass(frozen=True) class HaloWindow: - """Janela de leitura com halo, recortada nas bordas do domínio. + """Read window with halo, clipped at the domain edges. - global_slices : onde ler no array global (pode ser menor que - block+2*halo perto das bordas — não há padding em disco). - core_offset : deslocamento (linha, coluna) do início do "core" do - bloco dentro da janela lida, já que a janela pode começar mais - perto do core quando o halo foi recortado numa borda. + global_slices : where to read in the global array (may be smaller + than block+2*halo near the edges — there is no padding on disk). + core_offset : (row, col) offset of the start of the block's "core" + inside the window read, since the window may start closer to the + core when the halo was clipped at an edge. """ global_slices: tuple[slice, slice] core_offset: tuple[int, int] def halo_window(block: Block, shape: tuple[int, int], halo: int) -> HaloWindow: - """Calcula a janela com halo de um bloco, recortada nas bordas do domínio.""" + """Compute a block's halo window, clipped at the domain edges.""" height, width = shape r0 = max(0, block.r0 - halo) c0 = max(0, block.c0 - halo) @@ -80,17 +77,17 @@ def halo_window(block: Block, shape: tuple[int, int], halo: int) -> HaloWindow: class MemmapRasterWorkspace: """ - Armazena arrays 2D nomeados em disco (np.memmap), com decomposição - em blocos+halo, double-buffering e checkpoint de progresso. + Stores named 2-D arrays on disk (np.memmap), with block+halo + decomposition, double-buffering and progress checkpoints. - Genérico: não impõe nomes de variável nem domínio. Qualquer modelo - (dissmodel ou não) que opere sobre arrays nomeados 2D com regra de - vizinhança local pode usar este workspace. + Generic: it imposes no variable names or domain. Any model (dissmodel + or not) that operates on named 2-D arrays with a local neighbourhood + rule can use this workspace. - Uso típico - ---------- + Typical use + ----------- >>> ws = MemmapRasterWorkspace.create( - ... root=Path("/tmp/meu_workspace"), + ... root=Path("/tmp/my_workspace"), ... shape=(10000, 10000), ... arrays={"state": np.uint8}, ... block_h=512, block_w=512, halo=1, @@ -99,7 +96,7 @@ class MemmapRasterWorkspace: >>> for step in range(n_steps): ... for block in ws.blocks(): ... window = ws.read_block_with_halo(block) # dict[name, np.ndarray] - ... result = minha_regra(window) # dict[name, np.ndarray] + ... result = my_rule(window) # dict[name, np.ndarray] ... ws.write_block_core(block, result) ... ws.swap_buffers() ... ws.checkpoint(step) @@ -113,7 +110,7 @@ def __init__(self, root: Path) -> None: self.root = Path(root).resolve() metadata_path = self.root / self.METADATA if not metadata_path.is_file(): - raise FileNotFoundError(f"Workspace não inicializado: {self.root}") + raise FileNotFoundError(f"Workspace not initialized: {self.root}") self.metadata = json.loads(metadata_path.read_text(encoding="utf-8")) self.shape = tuple(int(v) for v in self.metadata["shape"]) self.block_h = int(self.metadata["block_h"]) @@ -151,37 +148,37 @@ def create( block_w: int, halo: int = 1, ) -> MemmapRasterWorkspace: - """Cria um workspace novo, com os dois slots do double-buffer. - - Os arquivos `.dat` nascem ESPARSOS: só as regiões efetivamente - escritas ocupam blocos no disco. Ler uma região nunca escrita - devolve zero — garantia do próprio sistema de arquivos POSIX, - idêntica ao que uma pré-escrita de zeros daria. Ou seja, a - semântica é a mesma de antes desta mudança; o que muda é só o - custo em disco (medido: um array de 4000x4000 float64 sai de - 122 MB reais para 0 MB até que algo seja escrito). - - Consequência prática, e limite desta economia: ela só aparece se - quem carrega DEIXAR blocos sem escrever. Um carregador que - preenche todo bloco — inclusive os vazios, com um sentinela como - NaN — torna o arquivo denso de novo. E deixar de escrever - significa que aquela região vale ZERO, não "ausente": para - domínios em que 0 é um valor válido (um código de classe, uma - elevação ao nível do mar) os dois casos ficam indistinguíveis. - - Distinguir "ausente" de "zero válido" precisaria de um canal a - mais, que este formato não tem — ver a nota "Esparsidade e custo - de disco" no README para o desenho proposto (índice de blocos - presentes + nodata declarado por array, o mesmo mecanismo que o - GeoTIFF usa com TileOffsets == 0). - - Nada disso afeta o custo de RAM, que é resolvido pelo memmap em - si: o kernel traz páginas sob demanda, então percorrer a grade - inteira nunca a materializa de uma vez. + """Create a new workspace, with both double-buffer slots. + + The `.dat` files start SPARSE: only regions actually written take + up blocks on disk. Reading a region never written returns zero — + guaranteed by the POSIX filesystem itself, identical to what + pre-writing zeros would give. The semantics are therefore + unchanged; only the disk cost differs (measured: a 4000x4000 + float64 array goes from 122 MB on disk to 0 MB until something + is written). + + Practical consequence, and the limit of this saving: it only + shows if the loader LEAVES blocks unwritten. A loader that fills + every block — including the empty ones, with a sentinel such as + NaN — makes the file dense again. And leaving a block unwritten + means that region is ZERO, not "missing": in domains where 0 is a + valid value (a class code, an elevation at sea level) the two + cases become indistinguishable. + + Telling "missing" apart from "valid zero" would need an extra + channel, which this format does not have — see "Sparsity and disk + cost" in the README for the proposed design (an index of present + blocks + nodata declared per array, the same mechanism GeoTIFF + uses with TileOffsets == 0). + + None of this affects RAM cost, which the memmap itself solves: + the kernel brings pages in on demand, so walking the whole grid + never materializes it at once. """ root = Path(root).resolve() if root.exists() and any(root.iterdir()): - raise FileExistsError(f"O diretório do workspace deve ser novo e vazio: {root}") + raise FileExistsError(f"The workspace directory must be new and empty: {root}") root.mkdir(parents=True, exist_ok=True) dtypes = {n: np.dtype(d) for n, d in arrays.items()} @@ -189,9 +186,9 @@ def create( for name, dtype in dtypes.items(): path = root / slot / f"{name}.dat" path.parent.mkdir(parents=True, exist_ok=True) - # mode="w+" dimensiona o arquivo sem tocar os bytes: ele fica - # esparso. NÃO pré-escrever zeros aqui — isso alocaria tudo - # fisicamente sem mudar nada do que se lê depois. + # mode="w+" sizes the file without touching its bytes: it stays + # sparse. Do NOT pre-write zeros here — that would allocate + # everything physically without changing anything read later. mm = np.memmap(path, dtype=dtype, mode="w+", shape=shape) mm.flush() @@ -205,29 +202,30 @@ def create( _write_json_atomic(root / cls.CHECKPOINT, {"step": 0, "read_slot": "a"}) return cls(root) - # ── acesso a blocos ────────────────────────────────────────────── + # ── block access ───────────────────────────────────────────────── def blocks(self) -> list[Block]: return make_blocks(self.shape[0], self.shape[1], self.block_h, self.block_w) def fill(self, name: str, array: np.ndarray, slot: str | None = None) -> None: - """Popula um array inteiro (uso único, ex.: estado inicial). Para - arrays grandes, prefira escrever bloco a bloco via write_block_core.""" + """Fill a whole array (one-off use, e.g. the initial state). For + large arrays, prefer writing block by block with write_block_core.""" slot = slot or self.checkpoint_data["read_slot"] mm = self._slots[slot][name] mm[:] = np.asarray(array, dtype=mm.dtype) mm.flush() def read_block_with_halo(self, block: Block, boundary_value=0) -> dict[str, np.ndarray]: - """Lê a janela com halo de todos os arrays para um bloco, do slot - de leitura atual. Preenche com boundary_value apenas o halo que - cai fora do domínio (bordas externas da grade) — o resto é lido - diretamente do disco, sem materializar a grade inteira. - - boundary_value aceita um escalar (mesmo valor para todos os - arrays) ou um dict {nome: valor} — necessário sempre que 0 não - for um sentinela seguro para algum array (ex.: um código de - classe válido no domínio). Ver engine.resolve_boundary_value.""" + """Read the halo window of every array for a block, from the + current read slot. Only the part of the halo that falls outside + the domain (the grid's outer edges) is filled with boundary_value + — the rest is read straight from disk, without materializing the + whole grid. + + boundary_value takes a scalar (same value for every array) or a + dict {name: value} — needed whenever 0 is not a safe sentinel for + some array (e.g. a valid class code in the domain). See + engine.resolve_boundary_value.""" read_slot = self._slots[self.checkpoint_data["read_slot"]] window = halo_window(block, self.shape, self.halo) full_h = block.r1 - block.r0 + 2 * self.halo @@ -240,8 +238,8 @@ def read_block_with_halo(self, block: Block, boundary_value=0) -> dict[str, np.n if raw.shape == (full_h, full_w): result[name] = raw.copy() else: - # bloco na borda: janela foi recortada, preenche o - # halo faltante com o valor de contorno daquele array. + # edge block: the window was clipped; fill the missing + # halo with that array's boundary value. padded = np.full((full_h, full_w), name_boundary, dtype=raw.dtype) dest_r0 = self.halo - window.core_offset[0] dest_c0 = self.halo - window.core_offset[1] @@ -250,51 +248,49 @@ def read_block_with_halo(self, block: Block, boundary_value=0) -> dict[str, np.n return result def write_block_core(self, block: Block, values: dict[str, np.ndarray]) -> None: - """Escreve a região core (sem halo) de um bloco no slot de - ESCRITA (o outro slot, não o de leitura) — preserva o - double-buffer: o passo atual nunca escreve no array que ainda - está sendo lido por outros blocos do mesmo passo.""" + """Write a block's core region (no halo) to the WRITE slot (the + other slot, not the read one) — keeps the double buffer intact: + the current step never writes to the array other blocks of the + same step are still reading.""" write_slot_name = "b" if self.checkpoint_data["read_slot"] == "a" else "a" write_slot = self._slots[write_slot_name] for name, core in values.items(): write_slot[name][block.core] = core def read_block_core(self, block: Block, name: str) -> np.ndarray: - """Lê só a região core (sem halo) de um array, do slot de - leitura atual. Usado para sincronização "_past" — não precisa - de vizinhança, só uma cópia direta.""" + """Read only the core region (no halo) of one array, from the + current read slot. Used for "_past" synchronization — no + neighbourhood needed, just a direct copy.""" read_slot = self._slots[self.checkpoint_data["read_slot"]] return np.asarray(read_slot[name][block.core]).copy() def write_block_to_read_slot(self, block: Block, name: str, values: np.ndarray) -> None: - """Escreve no MESMO slot que está sendo lido no momento (não no - slot de escrita do ping-pong). Uso exclusivo para popular - arrays "_past": eles devem estar disponíveis no slot de - leitura ANTES de execute() rodar naquele mesmo passo — ao - contrário dos arrays "correntes", que só ficam prontos no - próximo passo, após o swap.""" + """Write to the SAME slot currently being read (not to the + ping-pong write slot). Only for filling "_past" arrays: + they must be available in the read slot BEFORE execute() runs in + that same step — unlike the "current" arrays, which only become + ready in the next step, after the swap.""" read_slot = self._slots[self.checkpoint_data["read_slot"]] read_slot[name][block.core] = values def write_block_core_in_place(self, block: Block, values: dict[str, np.ndarray]) -> None: - """Escreve vários arrays no MESMO slot de leitura, imediatamente - visível a qualquer bloco processado depois (dentro da mesma - varredura). Sem ping-pong: propositalmente diferente de - write_block_core() (que escreve no slot oposto, só visível - após swap_buffers()). - - Uso: algoritmos de convergência iterativa (ver - sweep_until_convergence em convergence.py), onde não existe - "estado congelado do início do passo" — é refinamento - sucessivo do MESMO estado até um ponto fixo, e ver a atualização - do bloco vizinho processado momentos antes acelera a - convergência (iteração Gauss-Seidel, não Jacobi).""" + """Write several arrays to the SAME read slot, immediately + visible to any block processed later (within the same sweep). No + ping-pong: deliberately different from write_block_core() (which + writes to the opposite slot, visible only after swap_buffers()). + + Use: iterative convergence algorithms (see sweep_until_convergence + in convergence.py), where there is no "state frozen at the start + of the step" — it is successive refinement of the SAME state up + to a fixed point, and seeing the update of the neighbouring block + processed moments earlier speeds up convergence (Gauss-Seidel + iteration, not Jacobi).""" read_slot = self._slots[self.checkpoint_data["read_slot"]] for name, core in values.items(): read_slot[name][block.core] = core def swap_buffers(self) -> None: - """Troca leitura/escrita ao final de um passo de tempo completo.""" + """Swap read/write slots at the end of a complete time step.""" self.checkpoint_data["read_slot"] = ( "b" if self.checkpoint_data["read_slot"] == "a" else "a" ) @@ -304,9 +300,9 @@ def checkpoint(self, step: int) -> None: _write_json_atomic(self.root / self.CHECKPOINT, self.checkpoint_data) def snapshot(self, name: str) -> np.ndarray: - """Lê o array inteiro do slot de leitura atual (materializa em - RAM — usar só para inspeção/teste, nunca dentro do loop de - blocos).""" + """Read the whole array from the current read slot (materializes + it in RAM — for inspection/tests only, never inside the block + loop).""" return np.asarray(self._slots[self.checkpoint_data["read_slot"]][name]).copy() def flush(self) -> None: @@ -315,7 +311,7 @@ def flush(self) -> None: mm.flush() def as_backend(self, stride: int = 1, nodata_value: float | None = None): - """Devolve um adaptador WorkspaceRasterBackend para uso com RasterMap/dissmodel.""" + """Return a WorkspaceRasterBackend adapter for use with RasterMap/dissmodel.""" from .backend import WorkspaceRasterBackend return WorkspaceRasterBackend(self, stride=stride, nodata_value=nodata_value) diff --git a/src/haloexec/engine.py b/src/haloexec/engine.py index 912f823..d1a1bbf 100644 --- a/src/haloexec/engine.py +++ b/src/haloexec/engine.py @@ -1,16 +1,16 @@ """ -Primitivas genéricas de Decomposição de Domínio com zonas de Halo -(Ghost Cell Pattern) — sem dependência de dissmodel. +Generic Domain Decomposition primitives with Halo zones (Ghost Cell +Pattern) — no dependency on dissmodel. -Fundamentação teórica: Kjolstad & Snir (2010), "Ghost Cell Pattern", -ParaPLoP; aplicação em AC-LULC geoespacial: Xia et al. (2025), ISPRS -IJGI 14(3):109. Ver README.md. +Theoretical basis: Kjolstad & Snir (2010), "Ghost Cell Pattern", +ParaPLoP; application to geospatial CA-LULC: Xia et al. (2025), ISPRS +IJGI 14(3):109. See README.md. -Este módulo contém apenas a lógica de particionamento da grade -(Block, make_blocks). A execução em si — chamar a regra de transição -por bloco e reconciliar o resultado — é responsabilidade de quem -consome estas primitivas. A integração concreta com o dissmodel está -em `dissmodel_ca.py` (HaloChunkedRasterCellularAutomaton). +This module holds only the grid-partitioning logic (Block, make_blocks). +The execution itself — calling the transition rule per block and +reconciling the result — belongs to whoever consumes these primitives. +The concrete dissmodel integration is in `ram/cellular_automaton.py` +(HaloChunkedRasterCellularAutomaton). """ from __future__ import annotations @@ -20,7 +20,7 @@ @dataclass(frozen=True) class Block: - """Um sub-domínio retangular da grade global (sem halo).""" + """A rectangular sub-domain of the global grid (no halo).""" r0: int r1: int @@ -33,16 +33,16 @@ def shape(self) -> tuple[int, int]: @property def core(self) -> tuple[slice, slice]: - """Slices prontos pra indexar a região deste bloco num array - global (usado pela camada de disco, ex.: write_block_core).""" + """Slices ready to index this block's region in a global array + (used by the disk layer, e.g. write_block_core).""" return (slice(self.r0, self.r1), slice(self.c0, self.c1)) def make_blocks(height: int, width: int, block_h: int, block_w: int) -> list[Block]: - """Decompõe uma grade (height, width) em blocos regulares de tamanho - (block_h, block_w). Blocos na borda direita/inferior podem ser - menores (resíduo), conforme decomposição de domínio regular - (Xia et al. 2025, Seção 2.1).""" + """Split a (height, width) grid into regular blocks of size + (block_h, block_w). Blocks on the right/bottom edge may be smaller + (the remainder), as in regular domain decomposition + (Xia et al. 2025, Section 2.1).""" blocks = [] for r0 in range(0, height, block_h): r1 = min(r0 + block_h, height) @@ -53,24 +53,23 @@ def make_blocks(height: int, width: int, block_h: int, block_w: int) -> list[Blo def resolve_boundary_value(boundary_value, name: str) -> float: - """Resolve o valor de preenchimento do halo externo para um array - específico. Aceita um escalar (mesmo valor para todos os arrays) ou - um dict {nome: valor}. + """Resolve the outer-halo fill value for one array. Accepts a + scalar (same value for every array) or a dict {name: value}. - Se `name` terminar em "_past" e não tiver entrada própria no dict, - cai automaticamente para o valor do nome base (sem "_past") — assim - quem configura {"solo": -1} não precisa lembrar de duplicar para - "solo_past" também. Sem esse fallback, "_past" cairia - silenciosamente no default 0, reintroduzindo o mesmo problema que - este mecanismo existe para evitar. + If `name` ends in "_past" and has no entry of its own in the dict, + it falls back to the base name's value (without "_past") — so + whoever sets {"solo": -1} does not have to remember to repeat it for + "solo_past". Without this fallback, "_past" would silently get the + default 0, bringing back the very problem this mechanism exists to + prevent. - Importante: 0 não é um sentinela seguro para todo domínio — em - BR-MANGUE, por exemplo, `solo=0` é SOLO_CANAL_FLUVIAL, um código - de solo VÁLIDO (não "sem dado"). Usar 0 como boundary_value para - esse array cria fontes de migração fantasmas na borda externa do - domínio, divergindo do resultado monolítico. Prefira alinhar - boundary_value ao nodata real de cada array (ex.: TIFF_BANDS do - domínio), passando um dict em vez de um escalar único. + Important: 0 is not a safe sentinel for every domain — in BR-MANGUE, + for example, `solo=0` is SOLO_CANAL_FLUVIAL, a VALID soil code (not + "no data"). Using 0 as boundary_value for that array creates phantom + migration sources at the domain's outer edge, diverging from the + monolithic result. Prefer aligning boundary_value with each array's + real nodata (e.g. the domain's TIFF_BANDS), passing a dict instead of + a single scalar. """ if not isinstance(boundary_value, dict): return boundary_value diff --git a/src/haloexec/ram/cellular_automaton.py b/src/haloexec/ram/cellular_automaton.py index 52b82eb..ef64ac8 100644 --- a/src/haloexec/ram/cellular_automaton.py +++ b/src/haloexec/ram/cellular_automaton.py @@ -1,37 +1,35 @@ """ -Integração com dissmodel: HaloChunkedRasterCellularAutomaton. - -Estende dissmodel.geo.raster.cellular_automaton.RasterCellularAutomaton -para executar rule() em blocos com halo, em vez de sobre a grade -inteira de uma vez. - -Ponto central de design: mantém o MESMO contrato de rule() da classe -base (`rule(arrays: dict[str, np.ndarray]) -> dict[str, np.ndarray]`). -Isso significa que qualquer RasterCellularAutomaton já escrito para -dissmodel roda em blocos+halo apenas trocando a classe base — nenhuma -mudança na lógica da regra é necessária. É essa propriedade que torna -a migração futura do BR-MANGUE (`chunked_engine.py`) para este motor -uma troca estrutural, não uma reescrita. - -Como funciona -------------- -1. Tira um snapshot da grade global (equivalente a `self.backend.past`). -2. Preenche halo global (`np.pad`) em cada array. -3. Para cada bloco, monta um RasterBackend temporário só com a - sub-grade + halo daquele bloco, e troca `self.backend` para ele. - Isso é o que garante que chamadas internas da regra como - `self.backend.focal_sum_mask(...)` operem sobre a forma local - correta (RasterBackend.focal_sum_mask usa `self.shape` do backend - ativo) em vez da forma global. -4. Chama `self.rule(block_backend.snapshot())` — mesma assinatura de - sempre. -5. Recorta o halo do resultado (mantém só a região "core") e escreve - na posição correspondente da grade global nova. -6. Restaura `self.backend` para o backend global real e aplica as - atualizações. - -Fundamentação teórica: Kjolstad & Snir (2010), Ghost Cell Pattern -(ParaPLoP); Xia et al. (2025), ISPRS IJGI 14(3):109 — ver README.md. +dissmodel integration: HaloChunkedRasterCellularAutomaton. + +Extends dissmodel.geo.raster.cellular_automaton.RasterCellularAutomaton +to run rule() in blocks with a halo, instead of over the whole grid at +once. + +Central design point: it keeps the SAME rule() contract as the base +class (`rule(arrays: dict[str, np.ndarray]) -> dict[str, np.ndarray]`). +Any RasterCellularAutomaton already written for dissmodel therefore +runs in blocks+halo just by swapping the base class — no change to the +rule's logic is needed. This property is what makes moving an existing +model onto this engine a structural swap, not a rewrite. + +How it works +------------ +1. Take a snapshot of the global grid (equivalent to `self.backend.past`). +2. Pad each array with the global halo (`np.pad`). +3. For each block, build a temporary RasterBackend holding only that + block's sub-grid + halo, and point `self.backend` at it. This is + what makes calls inside the rule such as + `self.backend.focal_sum_mask(...)` operate on the correct local + shape (RasterBackend.focal_sum_mask uses the active backend's + `self.shape`) instead of the global one. +4. Call `self.rule(block_backend.snapshot())` — the usual signature. +5. Crop the halo from the result (keep only the "core" region) and + write it at the matching position of the new global grid. +6. Restore `self.backend` to the real global backend and apply the + updates. + +Theoretical basis: Kjolstad & Snir (2010), Ghost Cell Pattern +(ParaPLoP); Xia et al. (2025), ISPRS IJGI 14(3):109 — see README.md. """ from __future__ import annotations @@ -45,12 +43,12 @@ class HaloChunkedRasterCellularAutomaton(RasterCellularAutomaton): """ - RasterCellularAutomaton que processa a grade em blocos com halo, - em vez de de uma vez só. + RasterCellularAutomaton that processes the grid in blocks with a + halo, instead of all at once. - Uso: qualquer subclasse existente de RasterCellularAutomaton pode - trocar a herança para esta classe e ganhar decomposição de domínio - sem alterar `rule()`. + Usage: any existing RasterCellularAutomaton subclass can switch its + base class to this one and gain domain decomposition without + changing `rule()`. Examples -------- @@ -81,18 +79,18 @@ def setup( # type: ignore[override] Parameters ---------- backend : RasterBackend - Backend global compartilhado (mesma semântica da classe base). + Shared global backend (same semantics as the base class). block_h, block_w : int - Dimensões do bloco de processamento. + Processing block size. halo : int, optional - Raio da vizinhança da regra. Deve ser >= alcance máximo de - dependência espacial de um passo de tempo. Default 1 - (Moore/Von Neumann de vizinho imediato). + The rule's neighbourhood radius. Must be >= the maximum + spatial dependency reach of one time step. Default 1 + (immediate Moore/Von Neumann neighbours). boundary_value : float, optional - Valor de preenchimento do halo global nas bordas externas - da grade (fora do domínio simulado). Default 0. + Fill value of the global halo at the grid's outer edges + (outside the simulated domain). Default 0. state_attr : str, optional - Ver classe base. + See the base class. """ super().setup(backend=backend, state_attr=state_attr) self.block_h = block_h @@ -101,7 +99,7 @@ def setup( # type: ignore[override] self.boundary_value = boundary_value def _block_backend(self, padded: dict[str, np.ndarray], block: Block) -> RasterBackend: - """Monta um RasterBackend temporário com a sub-grade+halo do bloco.""" + """Build a temporary RasterBackend with the block's sub-grid + halo.""" h = self.halo block_shape = (block.r1 - block.r0 + 2 * h, block.c1 - block.c0 + 2 * h) temp = RasterBackend(shape=block_shape) @@ -112,11 +110,11 @@ def _block_backend(self, padded: dict[str, np.ndarray], block: Block) -> RasterB def execute(self) -> None: """ - Executa um passo de tempo processando a grade em blocos+halo. + Run one time step, processing the grid in blocks+halo. - Substitui o execute() da classe base (que chama rule() uma vez - sobre a grade inteira) por um loop de blocos, preservando o - mesmo contrato de rule() para quem escreve a regra. + Replaces the base class's execute() (which calls rule() once over + the whole grid) with a block loop, keeping the same rule() + contract for whoever writes the rule. """ real_backend = self.backend height, width = real_backend.shape @@ -127,7 +125,7 @@ def execute(self) -> None: name: np.pad(arr, h, mode="constant", constant_values=resolve_boundary_value(self.boundary_value, name)) for name, arr in global_snapshot.items() - if arr.ndim == 2 # arrays temporais (time, y, x) não são suportados aqui + if arr.ndim == 2 # temporal (time, y, x) arrays are not supported here } new_arrays: dict[str, np.ndarray] = { @@ -137,9 +135,9 @@ def execute(self) -> None: for block in make_blocks(height, width, self.block_h, self.block_w): block_backend = self._block_backend(padded, block) - # Troca temporária: garante que chamadas como - # self.backend.focal_sum_mask(...) dentro de rule() operem - # sobre a forma local do bloco, não a forma global. + # Temporary swap: makes calls such as + # self.backend.focal_sum_mask(...) inside rule() operate on + # the block's local shape, not the global one. self.backend = block_backend try: updates = self.rule(block_backend.snapshot()) diff --git a/src/haloexec/ram/sync_model.py b/src/haloexec/ram/sync_model.py index 84b16ac..bba995d 100644 --- a/src/haloexec/ram/sync_model.py +++ b/src/haloexec/ram/sync_model.py @@ -1,37 +1,38 @@ """ -Integração com dissmodel: HaloChunkedSyncRasterModel. - -Diferente de RasterCellularAutomaton (que expõe um hook rule(arrays) -dedicado), modelos baseados em SyncRasterModel/RasterModel — como o -FloodModel do BR-MANGUE — implementam a lógica científica diretamente -em execute(), lendo/escrevendo arrays nomeados no backend -(self.backend.arrays["alt"], self.backend.get("uso_past"), etc.) e -usando self.shape/self.shift/self.dirs herdados de RasterModel. - -Este módulo generaliza a mesma estratégia de chunking+halo para esse -padrão, via herança múltipla cooperativa (mixin): HaloChunkedSyncRasterModel -intercepta execute() e setup(), delega a lógica real para a subclasse -concreta via super().execute(), com self.backend/self.shape trocados -temporariamente para uma sub-grade local por bloco. - -Isso significa que NENHUMA linha de FloodModel (ou de qualquer outro -SyncRasterModel) precisa mudar — apenas a ordem de herança na -declaração da classe: +dissmodel integration: HaloChunkedSyncRasterModel. + +Unlike RasterCellularAutomaton (which exposes a dedicated rule(arrays) +hook), models based on SyncRasterModel/RasterModel — such as a flood +model — implement their scientific logic directly in execute(), reading +and writing named arrays on the backend (self.backend.arrays["alt"], +self.backend.get("uso_past"), etc.) and using self.shape/self.shift/ +self.dirs inherited from RasterModel. + +This module extends the same chunking+halo strategy to that pattern, +through cooperative multiple inheritance (a mixin): +HaloChunkedSyncRasterModel intercepts execute() and setup() and +delegates the real logic to the concrete subclass through +super().execute(), with self.backend/self.shape temporarily swapped for +a per-block local sub-grid. + +This means NOT A SINGLE LINE of the model (FloodModel or any other +SyncRasterModel) has to change — only the inheritance order in the +class declaration: class FloodModelHalo(HaloChunkedSyncRasterModel, FloodModel): pass -A ordem importa (MRO): o mixin deve vir primeiro, para que seu -execute()/setup() seja chamado antes, com super() delegando para a -lógica real de FloodModel.execute()/setup(). - -Limitação conhecida: pre_execute()/post_execute() de SyncRasterModel -(que fazem o snapshot "_past") NÃO são interceptados por este -mixin — continuam operando sobre o backend global real, fora do loop -de blocos. Isso é intencional: sincronizar "_past" é uma cópia simples -de array inteiro, sem dependência de vizinhança, então não precisa de -decomposição de domínio. O halo só é necessário dentro de execute(), -onde há leitura de vizinhos via self.shift. +The order matters (MRO): the mixin must come first, so its +execute()/setup() runs first, with super() delegating to the real +FloodModel.execute()/setup(). + +Known limitation: SyncRasterModel's pre_execute()/post_execute() (which +take the "_past" snapshot) are NOT intercepted by this mixin — +they keep operating on the real global backend, outside the block loop. +This is intentional: synchronizing "_past" is a plain whole-array copy +with no neighbourhood dependency, so it needs no domain decomposition. +The halo is only needed inside execute(), where neighbours are read +through self.shift. """ from __future__ import annotations @@ -44,35 +45,35 @@ class FloodModelHalo(HaloChunkedSyncRasterModel, FloodModel): class HaloChunkedSyncRasterModel: """ - Mixin que processa execute() de um RasterModel/SyncRasterModel em - blocos com halo, delegando a lógica científica para a próxima - classe na MRO via super(). + Mixin that runs a RasterModel/SyncRasterModel's execute() in blocks + with a halo, delegating the scientific logic to the next class in + the MRO through super(). - Parameters (setup, além dos que a subclasse concreta já aceita) - ------------------------------------------------------------------ + Parameters (setup, on top of those the concrete subclass accepts) + ----------------------------------------------------------------- block_h, block_w : int - Dimensões do bloco de processamento. + Processing block size. halo : int, optional - Raio da vizinhança usado pela regra (default 1). - - ATENÇÃO — halo NÃO é sempre igual ao raio nominal do shift - usado pela regra. Se a regra computa uma quantidade DERIVADA - de vizinhos (ex.: um "fluxo" que depende de quantos vizinhos - satisfazem uma condição) e depois lê essa quantidade derivada - DE UM VIZINHO (não do próprio valor bruto), a dependência real - é de 2 saltos, não 1 — halo=1 fica sutilmente errado perto de - fronteiras internas de bloco (não nas bordas do domínio, que - já são tratadas por boundary_value). Achado documentado em - tests/test_flood_model_halo_depth_regression.py: o FloodModel - do BR-MANGUE precisa de halo=2 por esse motivo exato - (fluxo_viz depende de viz_baixos do vizinho, que depende dos - vizinhos do vizinho). Ao adaptar uma regra nova, se os testes - de equivalência passarem com dado sintético simples mas - falharem em dado real/irregular, suspeite de dependência de - 2+ saltos antes de suspeitar de outra coisa. + Neighbourhood radius used by the rule (default 1). + + WARNING — halo is NOT always equal to the nominal shift radius + the rule uses. If the rule computes a quantity DERIVED from + neighbours (e.g. a "flow" that depends on how many neighbours + satisfy a condition) and then reads that derived quantity FROM A + NEIGHBOUR (not its own raw value), the real dependency is 2 hops, + not 1 — halo=1 is then subtly wrong near internal block + boundaries (not at the domain edges, which boundary_value already + handles). Documented in the README ("the correct halo depth is + the dependency chain's depth"): the BR-MANGUE FloodModel needs + halo=2 for exactly this reason (its neighbour flow depends on the + neighbour's count of lower cells, which depends on the + neighbour's neighbours). When adapting a new rule, if the + equivalence tests pass on simple synthetic data but fail on + real/irregular data, suspect a 2+ hop dependency before anything + else. boundary_value : float, optional - Valor de preenchimento do halo global nas bordas externas da - grade. Default 0. + Fill value of the global halo at the grid's outer edges. + Default 0. """ def setup(self, backend: RasterBackend, block_h: int, block_w: int, @@ -81,7 +82,7 @@ def setup(self, backend: RasterBackend, block_h: int, block_w: int, self.block_w = block_w self.halo = halo self.boundary_value = boundary_value - super().setup(backend=backend, **kwargs) # delega para a subclasse real + super().setup(backend=backend, **kwargs) # delegate to the real subclass def execute(self) -> None: real_backend = self.backend @@ -89,8 +90,8 @@ def execute(self) -> None: height, width = real_backend.shape h = self.halo - # Todos os arrays estáticos (2D) do backend global, sem distinção - # de nome — genérico o suficiente para qualquer modelo concreto. + # Every static (2-D) array of the global backend, whatever its + # name — generic enough for any concrete model. static_names = [n for n, a in real_backend.arrays.items() if a.ndim == 2] padded = { n: np.pad(real_backend.arrays[n], h, mode="constant", @@ -107,22 +108,22 @@ def execute(self) -> None: sub = arr[block.r0: block.r1 + 2 * h, block.c0: block.c1 + 2 * h] block_backend.set(name, sub) - # Troca temporária: garante que self.shape e self.backend - # (usados diretamente dentro de execute() da subclasse real, - # ex. `rows, cols = self.shape` no FloodModel) reflitam a - # forma local do bloco, não a forma global. + # Temporary swap: makes self.shape and self.backend (used + # directly inside the real subclass's execute(), e.g. + # `rows, cols = self.shape`) reflect the block's local shape, + # not the global one. self.backend = block_backend self.shape = block_backend.shape try: - super().execute() # lógica real (ex.: FloodModel.execute) + super().execute() # the real logic (e.g. FloodModel.execute) finally: self.backend = real_backend self.shape = real_shape - # Reconcilia: recorta o halo e escreve na grade global nova. - # Arrays "_past" são ignorados aqui — são geridos pelo - # synchronize() global em pre_execute()/post_execute(), não - # devem ser sobrescritos com fatias locais com halo. + # Reconcile: crop the halo and write into the new global grid. + # "_past" arrays are skipped here — they are managed by + # the global synchronize() in pre_execute()/post_execute() and + # must not be overwritten with local slices that carry a halo. for name, arr in block_backend.arrays.items(): if name.endswith("_past"): continue diff --git a/src/haloexec/visualization.py b/src/haloexec/visualization.py index 4d2c381..9b69e13 100644 --- a/src/haloexec/visualization.py +++ b/src/haloexec/visualization.py @@ -1,8 +1,9 @@ """ -Componentes de visualização e checkpoints para haloexec e dissmodel. +Visualization and checkpoint components for haloexec and dissmodel. -Fornece CheckpointRasterMap para desenhar e salvar quadros (PNG) em passos -ou anos específicos da simulação (evitando overhead em passos intermediários). +Provides CheckpointRasterMap, which draws and saves frames (PNG) only at +chosen simulation steps or years (avoiding the overhead on intermediate +steps). """ from __future__ import annotations @@ -20,20 +21,20 @@ if HAS_RASTERMAP: class CheckpointRasterMap(RasterMap): """ - Extensão do RasterMap que permite filtrar quais passos ou anos da simulação - serão desenhados e exportados para PNG. + RasterMap extension that filters which simulation steps or years are + drawn and exported to PNG. - Evita o processamento gráfico do matplotlib e a leitura de páginas de memmap - em passos intermediários, permitindo simulações de longa duração em grandes grades. + Skips matplotlib rendering and memmap page reads on intermediate steps, + making long simulations on large grids practical. - Parâmetros + Parameters ---------- save_steps : Iterable[int] | None - Lista ou conjunto de passos (ex.: [1, 5, 10, 20]) nos quais o quadro deve - ser renderizado e salvo. Se None, comporta-se como o RasterMap padrão - (obedecendo ao parâmetro `step` da classe base Model). + List or set of steps (e.g. [1, 5, 10, 20]) at which the frame is + rendered and saved. If None, behaves like the standard RasterMap + (following the `step` parameter of the base Model class). **kwargs - Todos os demais argumentos são repassados ao RasterMap (backend, band, + All other arguments are passed on to RasterMap (backend, band, color_map, cmap, save_frames, etc.). """ @@ -54,6 +55,6 @@ def execute(self) -> None: class CheckpointRasterMap: # type: ignore[no-redef] def __init__(self, *args: Any, **kwargs: Any) -> None: raise ImportError( - "RasterMap requer dissmodel instalado com extra viz: " + "RasterMap requires dissmodel installed with the viz extra: " "pip install 'dissmodel[viz]'" ) diff --git a/tests/test_convergence.py b/tests/test_convergence.py index ec92fbe..8adf758 100644 --- a/tests/test_convergence.py +++ b/tests/test_convergence.py @@ -1,12 +1,11 @@ """ -Prova de equivalência para sweep_until_convergence, usando exatamente -o caso de uso que motivou a primitiva: propagação de conectividade -(binary_propagation) por blocos+halo+varreduras, comparada a um -binary_propagation monolítico no domínio inteiro de uma vez. - -Isso é a prova que o protótipo original (chunked_engine.py::propagar_conectividade) -nunca teve — nenhum teste lá confirmava que a versão em blocos convergia -para o mesmo resultado exato que a versão monolítica. +Equivalence proof for sweep_until_convergence, using exactly the use +case that motivated the primitive: connectivity propagation +(binary_propagation) through blocks+halo+sweeps, compared with a +monolithic binary_propagation over the whole domain at once. + +It confirms that the block version converges to exactly the same result +as the monolithic one. """ import numpy as np @@ -19,13 +18,13 @@ def _connectivity_rule(window: dict[str, np.ndarray], halo: int = 1) -> dict[str, np.ndarray]: - """Mesma lógica do aluno: dilata 'conectado' através de 'permeavel', - dentro da janela com halo, e devolve só o núcleo.""" - connected = window["conectado"].astype(bool) - permeable = window["permeavel"].astype(bool) + """Dilate 'connected' through 'permeable' inside the halo window, + and return only the core.""" + connected = window["connected"].astype(bool) + permeable = window["permeable"].astype(bool) propagated = binary_propagation(connected, mask=permeable) core = propagated[halo:-halo, halo:-halo] - return {"conectado": core.astype(np.uint8)} + return {"connected": core.astype(np.uint8)} def _run_monolithic(seeds: np.ndarray, permeable: np.ndarray) -> np.ndarray: @@ -36,28 +35,28 @@ def _run_chunked(tmp_path, seeds: np.ndarray, permeable: np.ndarray, block_h: int, block_w: int, halo: int = 1) -> tuple[np.ndarray, dict]: ws = MemmapRasterWorkspace.create( root=tmp_path / "workspace", shape=seeds.shape, - arrays={"conectado": np.uint8, "permeavel": np.uint8}, + arrays={"connected": np.uint8, "permeable": np.uint8}, block_h=block_h, block_w=block_w, halo=halo, ) - ws.fill("conectado", seeds.astype(np.uint8)) - ws.fill("permeavel", permeable.astype(np.uint8)) + ws.fill("connected", seeds.astype(np.uint8)) + ws.fill("permeable", permeable.astype(np.uint8)) info = sweep_until_convergence( ws, lambda w: _connectivity_rule(w, halo), boundary_value=0, ) ws.flush() - return ws.snapshot("conectado"), info + return ws.snapshot("connected"), info def _labyrinth_scenario(height: int, width: int, seed: int) -> tuple[np.ndarray, np.ndarray]: - """Gera um labirinto de permeabilidade que força a conectividade a - serpentear por várias fronteiras de bloco antes de convergir — - testa de verdade a propagação através de múltiplos blocos, não só - vizinhança imediata de uma fonte central.""" + """Build a permeability labyrinth that forces connectivity to wind + across several block boundaries before converging — it really tests + propagation through many blocks, not just the immediate + neighbourhood of a central source.""" rng = np.random.default_rng(seed) - permeable = rng.random((height, width)) < 0.65 # a maioria é permeável + permeable = rng.random((height, width)) < 0.65 # most cells are permeable seeds = np.zeros((height, width), dtype=bool) - seeds[0, 0] = True # única fonte, no canto -- força propagação longa + seeds[0, 0] = True # single source, in the corner -- forces a long propagation permeable[0, 0] = True return seeds, permeable @@ -65,10 +64,10 @@ def _labyrinth_scenario(height: int, width: int, seed: int) -> tuple[np.ndarray, @pytest.mark.parametrize( "height, width, block_h, block_w, seed, label", [ - (40, 40, 10, 10, 42, "grade_divisivel_exatamente"), - (37, 53, 8, 12, 7, "grade_com_resto_blocos_irregulares"), - (30, 30, 6, 6, 123, "blocos_pequenos_muitas_fronteiras"), - (20, 20, 100, 100, 99, "bloco_maior_que_grade"), + (40, 40, 10, 10, 42, "grid_divides_exactly"), + (37, 53, 8, 12, 7, "grid_with_remainder_irregular_blocks"), + (30, 30, 6, 6, 123, "small_blocks_many_boundaries"), + (20, 20, 100, 100, 99, "block_larger_than_grid"), ], ) def test_sweep_until_convergence_equivalence(tmp_path, height, width, block_h, block_w, seed, label): @@ -78,7 +77,7 @@ def test_sweep_until_convergence_equivalence(tmp_path, height, width, block_h, b chunked, info = _run_chunked(tmp_path, seeds, permeable, block_h, block_w) n_diff = int(np.sum(golden != chunked)) - assert n_diff == 0, f"[{label}] {n_diff}/{height*width} células divergentes (info={info})" + assert n_diff == 0, f"[{label}] {n_diff}/{height*width} cells differ (info={info})" assert info["converged"] @@ -93,31 +92,31 @@ def test_sweep_until_convergence_stress_random_seeds(tmp_path, seed): def test_sweep_until_convergence_raises_if_never_converges(tmp_path): - """Regra que sempre 'muda' algo (nunca estabiliza) deve estourar - RuntimeError, não travar num loop silencioso.""" + """A rule that always 'changes' something (never settles) must raise + RuntimeError, not hang in a silent loop.""" ws = MemmapRasterWorkspace.create( root=tmp_path / "workspace", shape=(10, 10), - arrays={"contador": np.uint8}, block_h=5, block_w=5, halo=1, + arrays={"counter": np.uint8}, block_h=5, block_w=5, halo=1, ) - ws.fill("contador", np.zeros((10, 10), dtype=np.uint8)) + ws.fill("counter", np.zeros((10, 10), dtype=np.uint8)) - def regra_instavel(window): - # sempre incrementa -- nunca converge - core = window["contador"][1:-1, 1:-1] - return {"contador": (core + 1) % 250} + def unstable_rule(window): + # always increments -- never converges + core = window["counter"][1:-1, 1:-1] + return {"counter": (core + 1) % 250} - with pytest.raises(RuntimeError, match="não convergiu"): - sweep_until_convergence(ws, regra_instavel, max_sweeps=3) + with pytest.raises(RuntimeError, match="did not converge"): + sweep_until_convergence(ws, unstable_rule, max_sweeps=3) def test_sweep_until_convergence_reports_sweep_count(tmp_path): - """Uma única célula-fonte isolada (sem vizinho permeável) converge - na primeira varredura -- caso trivial, serve de sanity check do - contador de varreduras.""" + """A single isolated source cell (no permeable neighbour) converges + in the first sweep -- a trivial case, a sanity check of the sweep + counter.""" seeds = np.zeros((10, 10), dtype=bool) seeds[5, 5] = True permeable = np.zeros((10, 10), dtype=bool) - permeable[5, 5] = True # isolada -- nao tem pra onde propagar + permeable[5, 5] = True # isolated -- nowhere to propagate _, info = _run_chunked(tmp_path, seeds, permeable, block_h=5, block_w=5) assert info["sweeps"] == 1 diff --git a/tests/test_disk_backend_equivalence.py b/tests/test_disk_backend_equivalence.py index 5f3bb30..56e40b0 100644 --- a/tests/test_disk_backend_equivalence.py +++ b/tests/test_disk_backend_equivalence.py @@ -1,11 +1,11 @@ """ -Prova de equivalência do MemmapRasterWorkspace: execução em blocos+halo -lidos do DISCO (memmap, double-buffer, sem materializar a grade -inteira em memória) deve ser idêntica à execução monolítica em RAM. +Equivalence proof for MemmapRasterWorkspace: running in blocks+halo read +from DISK (memmap, double buffer, never materializing the whole grid in +memory) must match the monolithic in-RAM run exactly. -Este teste é deliberadamente independente de dissmodel — o workspace -em disco é genérico o suficiente para ser usado por qualquer framework -ou script solto, não só modelos dissmodel. +This test is deliberately independent of dissmodel — the on-disk +workspace is generic enough for any framework or standalone script, not +only dissmodel models. """ from pathlib import Path @@ -46,7 +46,7 @@ def _run_disk_backed(tmp_path: Path, initial: np.ndarray, generations: int, block_h=block_h, block_w=block_w, halo=halo, ) ws.fill("state", initial) - ws.fill("state", initial, slot="b") # ambos os slots começam iguais + ws.fill("state", initial, slot="b") # both slots start equal for step in range(generations): for block in ws.blocks(): diff --git a/tests/test_disk_backend_sparse.py b/tests/test_disk_backend_sparse.py index 249f869..f0bb0f7 100644 --- a/tests/test_disk_backend_sparse.py +++ b/tests/test_disk_backend_sparse.py @@ -1,16 +1,16 @@ """ -Alocação esparsa do MemmapRasterWorkspace. +Sparse allocation of MemmapRasterWorkspace. -`create()` dimensiona os `.dat` sem pré-escrever zeros, então eles nascem -esparsos. Estes testes fixam as DUAS metades do contrato: +`create()` sizes the `.dat` files without pre-writing zeros, so they +start sparse. These tests pin BOTH halves of the contract: - * a semântica não mudou — ler região nunca escrita continua devolvendo - zero, que é o que a pré-escrita dava; - * o custo em disco mudou — só o que é escrito ocupa blocos. + * the semantics are unchanged — reading a region never written still + returns zero, which is what pre-writing gave; + * the disk cost changed — only what is written takes up blocks. -O segundo depende do sistema de arquivos suportar arquivos esparsos (ext4, -xfs, btrfs suportam). Onde não suportar, o teste é pulado em vez de falhar: -a correção continua válida, só não rende ali. +The second depends on the filesystem supporting sparse files (ext4, xfs +and btrfs do). Where it does not, the test is skipped instead of +failing: the behaviour is still correct, it just saves nothing there. """ import numpy as np @@ -26,7 +26,7 @@ def _real_bytes(path): return path.stat().st_blocks * 512 -def _fs_suporta_esparso(tmp_path) -> bool: +def _fs_supports_sparse(tmp_path) -> bool: p = tmp_path / "probe.bin" mm = np.memmap(p, dtype="float64", mode="w+", shape=(1024, 1024)) mm.flush() @@ -41,64 +41,64 @@ def _ws(tmp_path, **kw): ) -# ── semântica: inalterada ──────────────────────────────────────────────────── +# ── semantics: unchanged ───────────────────────────────────────────────────── -def test_regiao_nunca_escrita_le_zero(tmp_path): - """Contrato antigo preservado: workspace novo lê zerado.""" +def test_region_never_written_reads_zero(tmp_path): + """Old contract kept: a new workspace reads as zeros.""" ws = _ws(tmp_path) for block in ws.blocks()[:5]: assert np.all(ws.read_block_core(block, "uso") == 0) -def test_halo_de_workspace_novo_le_zero(tmp_path): +def test_halo_of_new_workspace_reads_zero(tmp_path): ws = _ws(tmp_path) - janela = ws.read_block_with_halo(ws.blocks()[10], boundary_value=0)["uso"] - assert np.all(janela == 0) + window = ws.read_block_with_halo(ws.blocks()[10], boundary_value=0)["uso"] + assert np.all(window == 0) -def test_escrita_e_leitura_continuam_funcionando(tmp_path): +def test_write_and_read_still_work(tmp_path): ws = _ws(tmp_path) blk = ws.blocks()[3] - dados = np.full((blk.r1 - blk.r0, blk.c1 - blk.c0), 7.5) - ws.write_block_to_read_slot(blk, "uso", dados) + values = np.full((blk.r1 - blk.r0, blk.c1 - blk.c0), 7.5) + ws.write_block_to_read_slot(blk, "uso", values) ws.flush() - assert np.array_equal(ws.read_block_core(blk, "uso"), dados) + assert np.array_equal(ws.read_block_core(blk, "uso"), values) -def test_bloco_vizinho_de_um_escrito_continua_zero(tmp_path): - """Escrever um bloco não deve materializar valor em outro.""" +def test_neighbour_of_a_written_block_stays_zero(tmp_path): + """Writing one block must not materialize values in another.""" ws = _ws(tmp_path) - blocos = ws.blocks() + blocks = ws.blocks() ws.write_block_to_read_slot( - blocos[0], "uso", - np.ones((blocos[0].r1 - blocos[0].r0, blocos[0].c1 - blocos[0].c0)), + blocks[0], "uso", + np.ones((blocks[0].r1 - blocks[0].r0, blocks[0].c1 - blocks[0].c0)), ) ws.flush() - assert np.all(ws.read_block_core(blocos[1], "uso") == 0) + assert np.all(ws.read_block_core(blocks[1], "uso") == 0) -# ── custo em disco: esse é o ganho ─────────────────────────────────────────── +# ── disk cost: this is the gain ────────────────────────────────────────────── -def test_workspace_novo_nao_ocupa_disco(tmp_path): - if not _fs_suporta_esparso(tmp_path): - pytest.skip("sistema de arquivos não suporta arquivos esparsos") +def test_new_workspace_takes_no_disk(tmp_path): + if not _fs_supports_sparse(tmp_path): + pytest.skip("filesystem does not support sparse files") _ws(tmp_path) dat = tmp_path / "ws" / "a" / "uso.dat" - assert dat.stat().st_size == SHAPE[0] * SHAPE[1] * 8 # tamanho aparente cheio - assert _real_bytes(dat) < dat.stat().st_size // 10 # mas quase nada real + assert dat.stat().st_size == SHAPE[0] * SHAPE[1] * 8 # full apparent size + assert _real_bytes(dat) < dat.stat().st_size // 10 # but almost nothing real -def test_so_o_que_foi_escrito_ocupa_disco(tmp_path): - if not _fs_suporta_esparso(tmp_path): - pytest.skip("sistema de arquivos não suporta arquivos esparsos") +def test_only_what_was_written_takes_disk(tmp_path): + if not _fs_supports_sparse(tmp_path): + pytest.skip("filesystem does not support sparse files") ws = _ws(tmp_path) dat = tmp_path / "ws" / "a" / "uso.dat" - blocos = ws.blocks() - for blk in blocos[:4]: + blocks = ws.blocks() + for blk in blocks[:4]: ws.write_block_to_read_slot( blk, "uso", np.ones((blk.r1 - blk.r0, blk.c1 - blk.c0)) ) ws.flush() - ocupado = _real_bytes(dat) - assert ocupado > 0 # o que foi escrito conta - assert ocupado < dat.stat().st_size // 2 # o resto não + used = _real_bytes(dat) + assert used > 0 # what was written counts + assert used < dat.stat().st_size // 2 # the rest does not diff --git a/tests/test_equivalence.py b/tests/test_equivalence.py index 53b6be3..aa05a58 100644 --- a/tests/test_equivalence.py +++ b/tests/test_equivalence.py @@ -1,18 +1,18 @@ """ -Prova de equivalência matemática usando a maquinaria REAL do dissmodel -(Environment, RasterBackend, RasterCellularAutomaton) — não um harness -isolado. +Mathematical equivalence proof using the REAL dissmodel machinery +(Environment, RasterBackend, RasterCellularAutomaton) — not an isolated +harness. -A mesma regra de Game of Life é escrita uma única vez (via mixin) e -executada por duas classes base diferentes: +The same Game of Life rule is written once (through a mixin) and run by +two different base classes: - - GameOfLifeMono (dissmodel.RasterCellularAutomaton) — monolítica - - GameOfLifeHalo (haloexec.HaloChunkedRasterCellularAutomaton) — em blocos + - GameOfLifeMono (dissmodel.RasterCellularAutomaton) — monolithic + - GameOfLifeHalo (haloexec.HaloChunkedRasterCellularAutomaton) — in blocks -O resultado final deve ser IDÊNTICO célula a célula após N passos de -tempo, para qualquer decomposição de domínio válida. Isso prova que -trocar a classe base (a mudança real que o BR-MANGUE fará na migração) -não altera o resultado científico do modelo. +The final result must be IDENTICAL cell by cell after N time steps, for +any valid domain decomposition. This proves that swapping the base class +(the actual change a model makes to migrate) does not change the model's +scientific result. """ import numpy as np @@ -28,7 +28,7 @@ class GameOfLifeRuleMixin: - """Regra escrita uma única vez, reutilizada pelas duas classes base.""" + """Rule written once, reused by both base classes.""" def rule(self, arrays): state = arrays["state"] diff --git a/tests/test_gameoflife_from_geotiff.py b/tests/test_gameoflife_from_geotiff.py index 2178555..2c59bbd 100644 --- a/tests/test_gameoflife_from_geotiff.py +++ b/tests/test_gameoflife_from_geotiff.py @@ -1,16 +1,14 @@ """ -Prova de equivalência ponta a ponta: Game of Life carregado de um -GeoTIFF (simulando o resultado de um mosaico já materializado — ver -mosaic_io.py) direto para MemmapRasterWorkspace, executado em -blocos+halo, comparado a uma execução monolítica em RAM carregada do -MESMO arquivo. - -Fecha o ciclo que faltava: os testes anteriores validam (a) o motor de -blocos+halo com dado sintético em RAM/disco, e (b) o carregamento de -TIFF/mosaico isoladamente (round-trip, sem rodar nenhum modelo em -cima). Este teste roda um modelo de verdade a partir de um TIFF de -verdade, incluindo um caso "grande" (maior que os testes sintéticos -anteriores) para dar mais confiança de escala. +End-to-end equivalence proof: Game of Life loaded from a GeoTIFF +(standing in for an already materialized mosaic) straight into a +MemmapRasterWorkspace, run in blocks+halo, compared with a monolithic +in-RAM run loaded from the SAME file. + +Closes the loop: the other tests validate (a) the blocks+halo engine +with synthetic data in RAM/on disk, and (b) TIFF/mosaic loading on its +own (round trip, with no model running on top). This test runs a real +model from a real TIFF, including a "large" case (bigger than the +synthetic tests) for more confidence at scale. """ from pathlib import Path @@ -80,8 +78,8 @@ def _run_disk_from_tiff(path: Path, tmp_path: Path, generations: int, @pytest.mark.parametrize( "height, width, block_h, block_w, generations, seed, label", [ - (40, 40, 10, 10, 10, 42, "pequeno_grade_divisivel"), - (37, 53, 8, 12, 8, 7, "pequeno_com_resto"), + (40, 40, 10, 10, 10, 42, "small_grid_divides_exactly"), + (37, 53, 8, 12, 8, 7, "small_with_remainder"), ], ) def test_gameoflife_from_geotiff_equivalence(tmp_path, height, width, block_h, block_w, @@ -89,20 +87,19 @@ def test_gameoflife_from_geotiff_equivalence(tmp_path, height, width, block_h, b rng = np.random.default_rng(seed) initial = (rng.random((height, width)) < 0.35).astype(np.uint8) - tif_path = tmp_path / "mosaico.tif" + tif_path = tmp_path / "mosaic.tif" _write_geotiff(tif_path, initial) golden = _run_monolithic_from_tiff(tif_path, generations) disk = _run_disk_from_tiff(tif_path, tmp_path, generations, block_h, block_w) - assert np.array_equal(golden, disk), f"[{label}] divergência pós-TIFF" + assert np.array_equal(golden, disk), f"[{label}] mismatch after TIFF round trip" def test_gameoflife_from_large_geotiff(): - """Caso 'grande': TIFF de 2000x2000 (4 milhões de células), gerado - bloco a bloco (nunca materializado inteiro em RAM na geração), - lido bloco a bloco, rodado em blocos+halo — a mesma cadeia - mosaico->TIFF->disco->halo que seria usada com um mosaico real.""" + """'Large' case: a 2000x2000 TIFF (4 million cells), read block by + block and run in blocks+halo — the same mosaic->TIFF->disk->halo + chain that would be used with a real mosaic.""" import shutil import tempfile @@ -112,7 +109,7 @@ def test_gameoflife_from_large_geotiff(): rng = np.random.default_rng(99) initial = (rng.random((height, width)) < 0.35).astype(np.uint8) - tif_path = tmp_path / "mosaico_grande.tif" + tif_path = tmp_path / "large_mosaic.tif" _write_geotiff(tif_path, initial) generations = 3 @@ -120,6 +117,6 @@ def test_gameoflife_from_large_geotiff(): disk = _run_disk_from_tiff(tif_path, tmp_path, generations, block_h=256, block_w=256) n_diff = int(np.sum(golden != disk)) - assert n_diff == 0, f"{n_diff}/{height*width} células divergentes no caso grande" + assert n_diff == 0, f"{n_diff}/{height*width} cells differ in the large case" finally: shutil.rmtree(tmp_path, ignore_errors=True) diff --git a/tests/test_geomosaic_integration.py b/tests/test_geomosaic_integration.py index c6eb87a..b834fc3 100644 --- a/tests/test_geomosaic_integration.py +++ b/tests/test_geomosaic_integration.py @@ -1,12 +1,12 @@ """ -Teste de integração haloexec <-> geomosaic (pacotes separados, sem -dependência de runtime entre si). - -Prova que load_geotiff_into_workspace (haloexec) lê corretamente, -bloco a bloco, um mosaico produzido pelo geomosaic — inclusive para -blocos que cruzam a fronteira entre tiles. geomosaic é usado aqui só -como dependência de TESTE (extra "geomosaic"), nunca importado pelo -código de runtime do haloexec. +haloexec <-> geomosaic integration test (separate packages, no runtime +dependency between them). + +Proves that load_geotiff_into_workspace (haloexec) correctly reads, block +by block, a mosaic produced by geomosaic — including blocks that cross +the boundary between tiles. geomosaic is used here only as a TEST +dependency (the "geomosaic" extra), never imported by haloexec's runtime +code. """ import numpy as np @@ -48,8 +48,8 @@ def test_haloexec_reads_geomosaic_vrt_across_tile_boundary(tmp_path): contract = geomosaic.build_mosaic_contract(paths) vrt_path = geomosaic.write_vrt(contract, tmp_path / "mosaico.vrt") - # blocos de 6x6 cruzam a fronteira do mosaico (coluna/linha 10) - # várias vezes -- nenhum código especial de mosaico no haloexec. + # 6x6 blocks cross the mosaic boundary (column/row 10) several + # times -- no mosaic-specific code in haloexec. ws = MemmapRasterWorkspace.create( root=tmp_path / "workspace", shape=(20, 20), arrays={"estado": np.int16}, block_h=6, block_w=6, halo=1, diff --git a/tests/test_geotiff_io_equivalence.py b/tests/test_geotiff_io_equivalence.py index 20fc4b5..34ba29f 100644 --- a/tests/test_geotiff_io_equivalence.py +++ b/tests/test_geotiff_io_equivalence.py @@ -1,9 +1,9 @@ """ -Prova de equivalência: carregar um GeoTIFF bloco a bloco direto para -MemmapRasterWorkspace deve produzir exatamente os mesmos valores que -carregar o arquivo inteiro em memória (via rasterio puro, como -referência) — para cada array declarado, em toda a extensão do -domínio, incluindo blocos de borda com resto. +Equivalence proof: loading a GeoTIFF block by block straight into a +MemmapRasterWorkspace must produce exactly the same values as loading +the whole file into memory (with plain rasterio, as the reference) — for +every declared array, over the whole domain, including edge blocks with +a remainder. """ import numpy as np @@ -42,9 +42,9 @@ def _write_test_geotiff(path, height, width, seed): @pytest.mark.parametrize( "height, width, block_h, block_w, seed, label", [ - (40, 40, 10, 10, 42, "grade_divisivel_exatamente"), - (37, 53, 8, 12, 7, "grade_com_resto_blocos_irregulares"), - (25, 25, 100, 100, 99, "bloco_maior_que_grade"), + (40, 40, 10, 10, 42, "grid_divides_exactly"), + (37, 53, 8, 12, 7, "grid_with_remainder_irregular_blocks"), + (25, 25, 100, 100, 99, "block_larger_than_grid"), ], ) def test_load_geotiff_into_workspace_equivalence(tmp_path, height, width, block_h, block_w, seed, label): @@ -59,9 +59,9 @@ def test_load_geotiff_into_workspace_equivalence(tmp_path, height, width, block_ ) load_geotiff_into_workspace(ws, tif_path, BAND_SPEC) - assert np.array_equal(ws.snapshot("uso"), uso), f"[{label}] 'uso' divergente" - assert np.array_equal(ws.snapshot("solo"), solo), f"[{label}] 'solo' divergente" - assert np.allclose(ws.snapshot("alt"), alt, atol=1e-5), f"[{label}] 'alt' divergente" + assert np.array_equal(ws.snapshot("uso"), uso), f"[{label}] 'uso' differs" + assert np.array_equal(ws.snapshot("solo"), solo), f"[{label}] 'solo' differs" + assert np.allclose(ws.snapshot("alt"), alt, atol=1e-5), f"[{label}] 'alt' differs" def test_load_geotiff_into_workspace_shape_mismatch_raises(tmp_path): @@ -70,23 +70,23 @@ def test_load_geotiff_into_workspace_shape_mismatch_raises(tmp_path): ws = MemmapRasterWorkspace.create( root=tmp_path / "workspace", - shape=(30, 30), # shape errado de propósito + shape=(30, 30), # wrong shape on purpose arrays={"uso": np.int16, "alt": np.float32, "solo": np.int16}, block_h=10, block_w=10, halo=1, ) - with pytest.raises(ValueError, match="não bate"): + with pytest.raises(ValueError, match="does not match"): load_geotiff_into_workspace(ws, tif_path, BAND_SPEC) def test_load_geotiff_into_workspace_partial_bands(tmp_path): - """Só declarar 'uso' e 'alt' no workspace deve ignorar 'solo' sem erro.""" + """Declaring only 'uso' and 'alt' in the workspace must skip 'solo' without error.""" tif_path = tmp_path / "test.tif" uso, alt, _solo = _write_test_geotiff(tif_path, 15, 15, seed=5) ws = MemmapRasterWorkspace.create( root=tmp_path / "workspace", shape=(15, 15), - arrays={"uso": np.int16, "alt": np.float32}, # sem "solo" de propósito + arrays={"uso": np.int16, "alt": np.float32}, # no "solo" on purpose block_h=5, block_w=5, halo=1, ) load_geotiff_into_workspace(ws, tif_path, BAND_SPEC) diff --git a/tests/test_gol_patterns_example.py b/tests/test_gol_patterns_example.py index 8ff2103..eecd459 100644 --- a/tests/test_gol_patterns_example.py +++ b/tests/test_gol_patterns_example.py @@ -1,22 +1,19 @@ """ -Prova de equivalência para examples/gol_patterns/gol_patterns_haloexec.py: -padrões clássicos de Game of Life (glider, blinker, beacon, toad, -block, pulsar), posicionados deliberadamente sobre fronteiras de -bloco, devem produzir resultado IDÊNTICO entre a execução monolítica -(RasterCellularAutomaton puro) e em blocos+halo -(HaloChunkedRasterCellularAutomaton). - -Isso é o que de fato sustenta a alegação do exemplo ("um padrão que -atravessa a borda de um bloco e continua se comportando como deveria -é a evidência visual mais direta de que o halo sincroniza -corretamente") -- rodar sem erro não prova isso; só a comparação -célula a célula prova. - -Achado ao escrever este teste: as coordenadas originais do exemplo -tinham beacon (9,25, 4x4) caindo inteiro dentro da área do pulsar -(4,20, 13x13) -- sobreposição silenciosa, já que a função place() -sobrescreve por atribuição direta. Corrigido movendo beacon para a -coluna 33. +Equivalence proof for examples/gol/gol_patterns_haloexec.py: classic +Game of Life patterns (glider, blinker, beacon, toad, block, pulsar), +deliberately placed over block boundaries, must produce an IDENTICAL +result between the monolithic run (plain RasterCellularAutomaton) and +the blocks+halo run (HaloChunkedRasterCellularAutomaton). + +This is what actually backs the example's claim ("a pattern that crosses +a block edge and keeps behaving as it should is the most direct visual +evidence that the halo synchronizes correctly") -- running without +error does not prove it; only the cell-by-cell comparison does. + +Found while writing this test: the example's original coordinates put +the beacon (9,25, 4x4) entirely inside the pulsar's area (4,20, 13x13) +-- a silent overlap, since place() overwrites by direct assignment. +Fixed by moving the beacon to column 33. """ import numpy as np @@ -34,8 +31,8 @@ ROWS, COLS = 40, 40 GENERATIONS = 16 -# mesmas posições do exemplo, já corrigidas (ver docstring acima) -POSICOES = { +# same positions as the example, already fixed (see the docstring above) +POSITIONS = { "glider": (8, 8), "blinker": (20, 5), "beacon": (9, 33), @@ -51,25 +48,25 @@ def _place(grid: np.ndarray, pattern: list[list[int]], top: int, left: int) -> N grid[top:top + h, left:left + w] = arr -def _grade_inicial() -> np.ndarray: +def _initial_grid() -> np.ndarray: grid = np.zeros((ROWS, COLS), dtype=np.int8) - for nome, (top, left) in POSICOES.items(): - _place(grid, PATTERNS[nome], top, left) + for name, (top, left) in POSITIONS.items(): + _place(grid, PATTERNS[name], top, left) return grid -def test_nenhum_padrao_se_sobrepoe(): - """Confirma que as posições não colidem entre si -- se colidissem, - place() sobrescreveria um padrão sobre o outro silenciosamente, - sem erro nenhum (foi exatamente o bug encontrado com o beacon - original em (9,25), antes da correção).""" - ocupacao = np.zeros((ROWS, COLS), dtype=int) - for nome, (top, left) in POSICOES.items(): - arr = np.array(PATTERNS[nome]) +def test_no_pattern_overlaps(): + """Confirm the positions do not collide -- if they did, place() would + silently write one pattern over another, with no error at all + (exactly the bug found with the original beacon at (9,25), before + the fix).""" + occupancy = np.zeros((ROWS, COLS), dtype=int) + for name, (top, left) in POSITIONS.items(): + arr = np.array(PATTERNS[name]) h, w = arr.shape - regiao = ocupacao[top:top + h, left:left + w] - assert regiao.sum() == 0, f"{nome} em ({top},{left}) colide com outro padrão já posicionado" - ocupacao[top:top + h, left:left + w] += 1 + region = occupancy[top:top + h, left:left + w] + assert region.sum() == 0, f"{name} at ({top},{left}) collides with a pattern already placed" + occupancy[top:top + h, left:left + w] += 1 class _GameOfLifeRuleMixin: @@ -92,13 +89,13 @@ class _GoLHalo(_GameOfLifeRuleMixin, HaloChunkedRasterCellularAutomaton): @pytest.mark.parametrize( "block_h, block_w, halo, label", [ - (10, 10, 1, "bloco_10x10_mesmo_do_exemplo"), - (7, 13, 1, "bloco_irregular_nao_alinhado_aos_padroes"), - (5, 5, 2, "bloco_pequeno_halo_maior"), + (10, 10, 1, "block_10x10_same_as_example"), + (7, 13, 1, "irregular_block_not_aligned_to_patterns"), + (5, 5, 2, "small_block_larger_halo"), ], ) -def test_padroes_classicos_equivalencia(block_h, block_w, halo, label): - grid0 = _grade_inicial() +def test_classic_patterns_equivalence(block_h, block_w, halo, label): + grid0 = _initial_grid() backend_mono = raster_grid(rows=ROWS, cols=COLS, attrs={"state": grid0.copy()}) env_mono = Environment(start_time=0, end_time=GENERATIONS) @@ -110,7 +107,7 @@ def test_padroes_classicos_equivalencia(block_h, block_w, halo, label): env_halo = Environment(start_time=0, end_time=GENERATIONS) _GoLHalo(backend=backend_halo, block_h=block_h, block_w=block_w, halo=halo, state_attr="state") env_halo.run() - resultado = backend_halo.arrays["state"].copy() + result = backend_halo.arrays["state"].copy() - n_diff = int(np.sum(golden != resultado)) - assert n_diff == 0, f"[{label}] {n_diff}/{ROWS*COLS} células divergentes" + n_diff = int(np.sum(golden != result)) + assert n_diff == 0, f"[{label}] {n_diff}/{ROWS*COLS} cells differ" diff --git a/tests/test_visualization_adapter.py b/tests/test_visualization_adapter.py index 191901e..18d7016 100644 --- a/tests/test_visualization_adapter.py +++ b/tests/test_visualization_adapter.py @@ -1,5 +1,5 @@ """ -Testes para WorkspaceRasterBackend e CheckpointRasterMap. +Tests for WorkspaceRasterBackend and CheckpointRasterMap. """ from pathlib import Path @@ -78,7 +78,7 @@ def test_checkpoint_raster_map_filtering(tmp_path: Path, monkeypatch): env = Environment(start_time=1, end_time=5) _GoLHalo(workspace=ws, halo=1, boundary_value=0) - # Só salvar os passos 1 e 4 + # Save only steps 1 and 4 save_steps = [1, 4] CheckpointRasterMap( backend=adapter, @@ -125,7 +125,7 @@ def test_save_workspace_to_geotiff_roundtrip(tmp_path: Path): assert np.array_equal(ds.read(1), uso_data) assert np.allclose(ds.read(2), alt_data) - # Verifica também o método as_backend do workspace + # Also check the workspace's as_backend method backend = ws1.as_backend(stride=2) assert backend.shape == (30, 40) assert np.array_equal(backend.arrays["uso"], uso_data[::2, ::2]) diff --git a/tests/test_zarr_axis_order_regression.py b/tests/test_zarr_axis_order_regression.py index 5c28c73..f1212c7 100644 --- a/tests/test_zarr_axis_order_regression.py +++ b/tests/test_zarr_axis_order_regression.py @@ -1,23 +1,24 @@ """ -Regressão: load_zarr_into_workspace deve ler corretamente um Zarr cuja -ordem de eixos em disco NÃO é (y, x) -- cenário real, não hipotético. - -O disscube (VariableWriter) grava via -`da.to_dataset(name=var_name).to_zarr(...)`, e o próprio CubeClient.load() -do disscube faz `.transpose("y", "x")` defensivamente antes de usar -qualquer array carregado -- evidência de que a ordem de eixos em disco -NÃO é garantida como (y, x). - -O perigo: um array QUADRADO com eixos trocados (x, y) em vez de (y, x) -tem o MESMO shape nos dois casos -- a checagem de shape sozinha não -detecta a inversão. Sem correção, isso corrompe silenciosamente -linha/coluna, sem erro nenhum. - -Reproduzido aqui com xarray real, gravando exatamente como o -VariableWriter do disscube grava (to_dataset().to_zarr()), não com -zarr.create_array() direto -- para capturar o metadado real de -dimension_names que o xarray grava (campo nativo do Zarr v3), a -mesma informação que a correção em zarr_io.py usa para normalizar. +Regression: load_zarr_into_workspace must correctly read a Zarr store +whose on-disk axis order is NOT (y, x) -- a real scenario, not a +hypothetical one. + +disscube (VariableWriter) writes through +`da.to_dataset(name=var_name).to_zarr(...)`, and disscube's own +CubeClient.load() calls `.transpose("y", "x")` defensively before using +any loaded array -- evidence that the on-disk axis order is NOT +guaranteed to be (y, x). + +The danger: a SQUARE array with swapped (x, y) axes instead of (y, x) +has the SAME shape either way -- a shape check alone does not detect +the swap. Without the fix, rows/columns are silently corrupted, with no +error at all. + +Reproduced here with real xarray, writing exactly as disscube's +VariableWriter does (to_dataset().to_zarr()), not with a direct +zarr.create_array() -- to capture the real dimension_names metadata +xarray writes (the native Zarr v3 field), the same information the fix +in disk/io/zarr.py uses to normalize. """ import numpy as np @@ -30,10 +31,10 @@ def _write_like_disscube(path, data_yx: np.ndarray, dims: tuple[str, ...], var_name: str): - """Grava exatamente como VariableWriter do disscube: + """Write exactly as disscube's VariableWriter does: da.to_dataset(name=...).to_zarr(..., mode="w", consolidated=False). - `dims` controla a ordem de eixos gravada em disco -- ("y","x") é o - caso "correto"/esperado, ("x","y") é o caso perigoso real.""" + `dims` sets the axis order written to disk -- ("y","x") is the + "correct"/expected case, ("x","y") the real dangerous one.""" if dims == ("y", "x"): raw = data_yx elif dims == ("x", "y"): @@ -45,10 +46,10 @@ def _write_like_disscube(path, data_yx: np.ndarray, dims: tuple[str, ...], var_n def test_load_zarr_handles_yx_axis_order(tmp_path): - """Caso 'correto' (y, x) -- deve continuar funcionando como sempre.""" + """The 'correct' (y, x) case -- must keep working as always.""" n = 6 data = np.arange(n * n).reshape(n, n).astype("int16") - store = tmp_path / "correto.zarr" + store = tmp_path / "correct.zarr" _write_like_disscube(store, data, ("y", "x"), "uso") ws = MemmapRasterWorkspace.create( @@ -61,21 +62,21 @@ def test_load_zarr_handles_yx_axis_order(tmp_path): def test_load_zarr_handles_xy_axis_order_square_array(tmp_path): - """Caso PERIGOSO: array QUADRADO gravado com eixos (x, y) -- - mesmo shape do caso correto, mas dado fisicamente transposto em - disco. Sem a correção de dimension_names, isso passaria a checagem - de shape e corromperia silenciosamente linha/coluna.""" + """DANGEROUS case: a SQUARE array written with (x, y) axes -- same + shape as the correct case, but the data is physically transposed on + disk. Without the dimension_names fix, this would pass the shape + check and silently corrupt rows/columns.""" n = 6 - # valores distintos por linha E coluna, para que uma transposição - # incorreta produza um array MENSURAVELMENTE diferente do original + # distinct values per row AND column, so a wrong transposition + # produces an array MEASURABLY different from the original data = np.arange(n * n).reshape(n, n).astype("int16") - store = tmp_path / "perigoso.zarr" + store = tmp_path / "swapped.zarr" _write_like_disscube(store, data, ("x", "y"), "uso") - # confirma que o shape em disco é IGUAL ao esperado (é exatamente - # isso que torna o bug silencioso sem a correção de eixo) + # confirm the on-disk shape EQUALS the expected one (exactly what + # makes the bug silent without the axis fix) root = zarr.open(str(store), mode="r") - assert root["uso"].shape == data.shape, "pré-condição do teste: shapes devem coincidir" + assert root["uso"].shape == data.shape, "test precondition: shapes must match" ws = MemmapRasterWorkspace.create( root=tmp_path / "workspace", shape=(n, n), @@ -83,24 +84,24 @@ def test_load_zarr_handles_xy_axis_order_square_array(tmp_path): ) load_zarr_into_workspace(ws, str(store), variable_map={"uso": "uso"}) - resultado = ws.snapshot("uso") - assert np.array_equal(resultado, data), ( - "load_zarr_into_workspace leu o array com linha/coluna trocadas -- " - "regressão do bug de ordem de eixos (x,y) vs (y,x)" + result = ws.snapshot("uso") + assert np.array_equal(result, data), ( + "load_zarr_into_workspace read the array with rows/columns swapped -- " + "regression of the (x,y) vs (y,x) axis-order bug" ) def test_load_zarr_handles_txy_axis_order_temporal(tmp_path): - """Variável temporal (3D) com ordem de eixos não-canônica: (x, y, time) - em vez de (time, y, x).""" + """Temporal (3-D) variable with a non-canonical axis order: (x, y, time) + instead of (time, y, x).""" n = 5 n_time = 3 - # (time, y, x) -- valores originais - serie = np.arange(n_time * n * n).reshape(n_time, n, n).astype("int16") - # grava fisicamente em ordem (x, y, time) - raw_disco = np.transpose(serie, (2, 1, 0)) # de (time,y,x) para (x,y,time) + # (time, y, x) -- original values + series = np.arange(n_time * n * n).reshape(n_time, n, n).astype("int16") + # physically written in (x, y, time) order + raw_on_disk = np.transpose(series, (2, 1, 0)) # from (time,y,x) to (x,y,time) - da = xr.DataArray(raw_disco, dims=("x", "y", "time"), name="mangue") + da = xr.DataArray(raw_on_disk, dims=("x", "y", "time"), name="mangue") store = tmp_path / "temporal.zarr" da.to_dataset(name="mangue").to_zarr(str(store), mode="w", consolidated=False) @@ -110,15 +111,15 @@ def test_load_zarr_handles_txy_axis_order_temporal(tmp_path): ) load_zarr_into_workspace(ws, str(store), variable_map={"mangue": "mangue"}, time_index=1) - assert np.array_equal(ws.snapshot("mangue"), serie[1]) + assert np.array_equal(ws.snapshot("mangue"), series[1]) def test_load_zarr_handles_xy_axis_order_zarr_v2_format(tmp_path): - """Store no FORMATO Zarr v2 (ex.: gravado por xarray/disscube mais - antigos, ou com zarr_format=2): não existe metadata.dimension_names, - o xarray guarda os nomes no atributo `_ARRAY_DIMENSIONS`. Mesmo - lendo com zarr-python 3, sem o fallback para esse atributo o array - quadrado (x, y) era carregado TRANSPOSTO, silenciosamente.""" + """A store in Zarr v2 FORMAT (e.g. written by older xarray/disscube, + or with zarr_format=2): there is no metadata.dimension_names; xarray + keeps the names in the `_ARRAY_DIMENSIONS` attribute. Even when read + with zarr-python 3, without the fallback to that attribute the square + (x, y) array was loaded TRANSPOSED, silently.""" n = 6 data = np.arange(n * n).reshape(n, n).astype("int16") store = tmp_path / "v2.zarr" @@ -132,6 +133,6 @@ def test_load_zarr_handles_xy_axis_order_zarr_v2_format(tmp_path): load_zarr_into_workspace(ws, str(store), variable_map={"uso": "uso"}) assert np.array_equal(ws.snapshot("uso"), data), ( - "Zarr formato v2 com eixos (x, y) carregado transposto -- " - "fallback para _ARRAY_DIMENSIONS ausente" + "Zarr v2-format store with (x, y) axes loaded transposed -- " + "missing fallback to _ARRAY_DIMENSIONS" ) diff --git a/tests/test_zarr_io.py b/tests/test_zarr_io.py index e9f5afd..a646c35 100644 --- a/tests/test_zarr_io.py +++ b/tests/test_zarr_io.py @@ -1,7 +1,7 @@ """ -Testes de zarr_io.py: carregamento de Zarr (grupo multi-variável, -array único, e com dimensão temporal — padrão disscube) direto para -MemmapRasterWorkspace, bloco a bloco. +Tests for disk/io/zarr.py: loading Zarr (multi-variable group, single +array, and with a time dimension — the disscube layout) straight into a +MemmapRasterWorkspace, block by block. """ import numpy as np @@ -13,13 +13,13 @@ def test_load_zarr_group_multi_variable(tmp_path): - """Simula o layout do disscube: um grupo com várias variáveis - (ex.: 'uso', 'alt'), cada uma acessada por nome.""" + """Mimic the disscube layout: a group with several variables + (e.g. 'uso', 'alt'), each accessed by name.""" rng = np.random.default_rng(1) uso = rng.integers(1, 9, size=(20, 20)).astype("int16") alt = rng.uniform(-2.0, 8.0, size=(20, 20)).astype("float32") - store_path = str(tmp_path / "grupo.zarr") + store_path = str(tmp_path / "group.zarr") root = zarr.open_group(store_path, mode="w") root.create_array("uso", shape=(20, 20), dtype="int16") root["uso"][:] = uso @@ -37,52 +37,52 @@ def test_load_zarr_group_multi_variable(tmp_path): def test_load_zarr_group_with_variable_map(tmp_path): - """Nome do array no workspace difere do nome da variável no zarr.""" + """The workspace array name differs from the variable name in the store.""" rng = np.random.default_rng(2) - dado = rng.integers(0, 100, size=(15, 15)).astype("int32") + data = rng.integers(0, 100, size=(15, 15)).astype("int32") - store_path = str(tmp_path / "grupo.zarr") + store_path = str(tmp_path / "group.zarr") root = zarr.open_group(store_path, mode="w") - root.create_array("dist_sedes", shape=(15, 15), dtype="int32") - root["dist_sedes"][:] = dado + root.create_array("dist_to_towns", shape=(15, 15), dtype="int32") + root["dist_to_towns"][:] = data ws = MemmapRasterWorkspace.create( root=tmp_path / "workspace", shape=(15, 15), - arrays={"distancia": np.int32}, block_h=5, block_w=5, halo=1, + arrays={"distance": np.int32}, block_h=5, block_w=5, halo=1, ) - load_zarr_into_workspace(ws, store_path, variable_map={"distancia": "dist_sedes"}) + load_zarr_into_workspace(ws, store_path, variable_map={"distance": "dist_to_towns"}) - assert np.array_equal(ws.snapshot("distancia"), dado) + assert np.array_equal(ws.snapshot("distance"), data) def test_load_zarr_single_array(tmp_path): - """Store é um único array (sem grupo).""" + """The store is a single array (no group).""" rng = np.random.default_rng(3) - dado = rng.integers(0, 10, size=(12, 12)).astype("uint8") + data = rng.integers(0, 10, size=(12, 12)).astype("uint8") - store_path = str(tmp_path / "unico.zarr") + store_path = str(tmp_path / "single.zarr") za = zarr.open_array(store_path, mode="w", shape=(12, 12), dtype="uint8") - za[:] = dado + za[:] = data ws = MemmapRasterWorkspace.create( root=tmp_path / "workspace", shape=(12, 12), - arrays={"estado": np.uint8}, block_h=4, block_w=4, halo=1, + arrays={"state": np.uint8}, block_h=4, block_w=4, halo=1, ) - load_zarr_into_workspace(ws, store_path, variable_map={"estado": None}) + load_zarr_into_workspace(ws, store_path, variable_map={"state": None}) - assert np.array_equal(ws.snapshot("estado"), dado) + assert np.array_equal(ws.snapshot("state"), data) def test_load_zarr_temporal_variable(tmp_path): - """Variável 3D (time, y, x) — padrão do 'Temporal Backend' do - disscube para produtos derivados com janela de validade.""" + """3-D variable (time, y, x) — the layout of disscube's 'Temporal + Backend' for derived products with a validity window.""" rng = np.random.default_rng(4) - serie = rng.integers(0, 5, size=(3, 10, 10)).astype("int16") # 3 anos + series = rng.integers(0, 5, size=(3, 10, 10)).astype("int16") # 3 years store_path = str(tmp_path / "temporal.zarr") root = zarr.open_group(store_path, mode="w") root.create_array("mangue", shape=(3, 10, 10), dtype="int16") - root["mangue"][:] = serie + root["mangue"][:] = series ws = MemmapRasterWorkspace.create( root=tmp_path / "workspace", shape=(10, 10), @@ -90,7 +90,7 @@ def test_load_zarr_temporal_variable(tmp_path): ) load_zarr_into_workspace(ws, store_path, time_index=1) - assert np.array_equal(ws.snapshot("mangue"), serie[1]) + assert np.array_equal(ws.snapshot("mangue"), series[1]) def test_load_zarr_temporal_without_time_index_raises(tmp_path): @@ -108,14 +108,14 @@ def test_load_zarr_temporal_without_time_index_raises(tmp_path): def test_load_zarr_shape_mismatch_raises(tmp_path): - store_path = str(tmp_path / "grupo.zarr") + store_path = str(tmp_path / "group.zarr") root = zarr.open_group(store_path, mode="w") root.create_array("uso", shape=(30, 30), dtype="int16") root["uso"][:] = np.zeros((30, 30), dtype="int16") ws = MemmapRasterWorkspace.create( - root=tmp_path / "workspace", shape=(20, 20), # shape errado de proposito + root=tmp_path / "workspace", shape=(20, 20), # wrong shape on purpose arrays={"uso": np.int16}, block_h=5, block_w=5, halo=1, ) - with pytest.raises(ValueError, match="não bate"): + with pytest.raises(ValueError, match="does not match"): load_zarr_into_workspace(ws, store_path) diff --git a/tests/test_zarr_tiles_integration.py b/tests/test_zarr_tiles_integration.py index c635468..e0da270 100644 --- a/tests/test_zarr_tiles_integration.py +++ b/tests/test_zarr_tiles_integration.py @@ -1,21 +1,22 @@ """ -Teste de integração haloexec <-> disscube (pacotes separados, sem -dependência de runtime entre si). - -Prova que `load_zarr_tiles_into_workspace` monta corretamente um mosaico -de N stores Zarr posicionados lado a lado — inclusive para blocos cujo -halo CRUZA a fronteira entre dois arquivos, e para buracos da malha. - -É o análogo, para Zarr, do que `test_geomosaic_integration.py` prova para -GeoTIFF/VRT. A diferença importa: no caso VRT o GDAL costura antes do -haloexec entrar em cena, então um erro de posicionamento apareceria já na -leitura do VRT. Aqui a costura é feita por este módulo, tile a tile — um -offset errado, ou um tile faltando, produz dado errado numa fronteira e -em mais lugar nenhum. Por isso os testes miram exatamente as fronteiras. - -O formato do layout (chaves e semântica) é o contrato com -`CubeClient.tile_layout()` do disscube; `test_layout_shape_matches_disscube` -o fixa dos dois lados quando o disscube está instalado. +haloexec <-> disscube integration test (separate packages, no runtime +dependency between them). + +Proves that `load_zarr_tiles_into_workspace` correctly assembles a +mosaic of N Zarr stores placed side by side — including blocks whose +halo CROSSES the boundary between two files, and holes in the tiling. + +It is the Zarr counterpart of what `test_geomosaic_integration.py` +proves for GeoTIFF/VRT. The difference matters: with a VRT, GDAL does +the stitching before haloexec comes in, so a positioning error would +already show when reading the VRT. Here the stitching is done by this +module, tile by tile — a wrong offset, or a missing tile, produces wrong +data at one boundary and nowhere else. That is why the tests aim +precisely at the boundaries. + +The layout format (keys and semantics) is the contract with disscube's +`CubeClient.tile_layout()`; `test_layout_shape_matches_disscube` pins it +from both sides when disscube is installed. """ import numpy as np @@ -29,12 +30,12 @@ load_zarr_tiles_into_workspace, ) -T = 8 # lado de cada tile — pequeno para as fronteiras ficarem inspecionáveis +T = 8 # side of each tile — small, so the boundaries stay inspectable def _write_tile(path, data): - """Grava um array 2D como store Zarr com dimension_names (y, x), - do mesmo jeito que xarray.to_zarr grava — que é como o disscube grava.""" + """Write a 2-D array as a Zarr store with dimension_names (y, x), + the same way xarray.to_zarr writes it — which is how disscube writes.""" root = zarr.open_group(str(path), mode="w") arr = root.create_array( "v", shape=data.shape, dtype=str(data.dtype), dimension_names=("y", "x") @@ -50,12 +51,12 @@ def _ws(tmp_path, shape, block=4, halo=1, dtype="float64"): ) -def _quadrantes(tmp_path, valores): - """Quatro tiles TxT num workspace 2Tx2T. `valores` dá o valor constante - de cada quadrante: (NO, NE, SO, SE). None = tile ausente (buraco).""" +def _quadrants(tmp_path, values): + """Four TxT tiles in a 2Tx2T workspace. `values` gives the constant + value of each quadrant: (NW, NE, SW, SE). None = missing tile (hole).""" pos = [(0, 0), (0, T), (T, 0), (T, T)] tiles = [] - for i, ((r, c), val) in enumerate(zip(pos, valores)): + for i, ((r, c), val) in enumerate(zip(pos, values)): if val is None: continue url = _write_tile(tmp_path / f"t{i}.zarr", np.full((T, T), val, dtype="float64")) @@ -66,123 +67,124 @@ def _quadrantes(tmp_path, valores): return tiles -# ── montagem básica ────────────────────────────────────────────────────────── +# ── basic assembly ─────────────────────────────────────────────────────────── def test_four_tiles_land_in_their_own_quadrants(tmp_path): ws = _ws(tmp_path, (2 * T, 2 * T)) - load_zarr_tiles_into_workspace(ws, _quadrantes(tmp_path, (1.0, 2.0, 3.0, 4.0))) + load_zarr_tiles_into_workspace(ws, _quadrants(tmp_path, (1.0, 2.0, 3.0, 4.0))) - lido = MemmapRasterWorkspace(tmp_path / "ws") + reopened = MemmapRasterWorkspace(tmp_path / "ws") full = np.empty((2 * T, 2 * T)) - for b in lido.blocks(): - full[b.r0:b.r1, b.c0:b.c1] = lido.read_block_core(b, "v") + for b in reopened.blocks(): + full[b.r0:b.r1, b.c0:b.c1] = reopened.read_block_core(b, "v") - assert np.all(full[:T, :T] == 1.0), "quadrante NO" - assert np.all(full[:T, T:] == 2.0), "quadrante NE" - assert np.all(full[T:, :T] == 3.0), "quadrante SO" - assert np.all(full[T:, T:] == 4.0), "quadrante SE" + assert np.all(full[:T, :T] == 1.0), "NW quadrant" + assert np.all(full[:T, T:] == 2.0), "NE quadrant" + assert np.all(full[T:, :T] == 3.0), "SW quadrant" + assert np.all(full[T:, T:] == 4.0), "SE quadrant" def test_row_and_column_offsets_are_not_swapped(tmp_path): - """Trocar row_off por col_off transporia o mosaico — e com tiles - quadrados o shape continuaria certo, então só o conteúdo denuncia.""" + """Swapping row_off and col_off would transpose the mosaic — and with + square tiles the shape would still be right, so only the content + gives it away.""" ws = _ws(tmp_path, (2 * T, 2 * T)) - load_zarr_tiles_into_workspace(ws, _quadrantes(tmp_path, (1.0, 2.0, 3.0, 4.0))) - lido = MemmapRasterWorkspace(tmp_path / "ws") - ne = lido.read_block_core(Block(r0=0, r1=4, c0=T, c1=T + 4), "v") - so = lido.read_block_core(Block(r0=T, r1=T + 4, c0=0, c1=4), "v") + load_zarr_tiles_into_workspace(ws, _quadrants(tmp_path, (1.0, 2.0, 3.0, 4.0))) + reopened = MemmapRasterWorkspace(tmp_path / "ws") + ne = reopened.read_block_core(Block(r0=0, r1=4, c0=T, c1=T + 4), "v") + so = reopened.read_block_core(Block(r0=T, r1=T + 4, c0=0, c1=4), "v") assert np.all(ne == 2.0) and np.all(so == 3.0) -# ── fronteiras: o que este módulo pode quebrar sozinho ─────────────────────── +# ── boundaries: what this module can break on its own ─────────────────────── def test_halo_crosses_boundary_between_two_zarr_files(tmp_path): - """A janela com halo tem de trazer o valor do tile VIZINHO — que veio - de outro arquivo — e não o do próprio tile nem nodata.""" + """The halo window must bring the value of the NEIGHBOURING tile — + which came from another file — not the tile's own value nor nodata.""" ws = _ws(tmp_path, (2 * T, 2 * T), block=4, halo=1) - load_zarr_tiles_into_workspace(ws, _quadrantes(tmp_path, (1.0, 2.0, 3.0, 4.0))) - lido = MemmapRasterWorkspace(tmp_path / "ws") + load_zarr_tiles_into_workspace(ws, _quadrants(tmp_path, (1.0, 2.0, 3.0, 4.0))) + reopened = MemmapRasterWorkspace(tmp_path / "ws") - # bloco encostado na fronteira vertical: halo à direita cai no tile NE + # block touching the vertical boundary: its right halo falls in the NE tile blk = Block(r0=0, r1=4, c0=T - 4, c1=T) - jan = lido.read_block_with_halo(blk, boundary_value=np.nan)["v"] - assert np.all(jan[1:-1, 1:-1] == 1.0), "núcleo é do tile NO" - assert np.all(jan[1:-1, -1] == 2.0), "halo direito tem de vir do tile NE" + win = reopened.read_block_with_halo(blk, boundary_value=np.nan)["v"] + assert np.all(win[1:-1, 1:-1] == 1.0), "the core belongs to the NW tile" + assert np.all(win[1:-1, -1] == 2.0), "the right halo must come from the NE tile" def test_halo_crosses_horizontal_boundary(tmp_path): ws = _ws(tmp_path, (2 * T, 2 * T), block=4, halo=1) - load_zarr_tiles_into_workspace(ws, _quadrantes(tmp_path, (1.0, 2.0, 3.0, 4.0))) - lido = MemmapRasterWorkspace(tmp_path / "ws") + load_zarr_tiles_into_workspace(ws, _quadrants(tmp_path, (1.0, 2.0, 3.0, 4.0))) + reopened = MemmapRasterWorkspace(tmp_path / "ws") blk = Block(r0=T - 4, r1=T, c0=0, c1=4) - jan = lido.read_block_with_halo(blk, boundary_value=np.nan)["v"] - assert np.all(jan[1:-1, 1:-1] == 1.0) - assert np.all(jan[-1, 1:-1] == 3.0), "halo inferior tem de vir do tile SO" + win = reopened.read_block_with_halo(blk, boundary_value=np.nan)["v"] + assert np.all(win[1:-1, 1:-1] == 1.0) + assert np.all(win[-1, 1:-1] == 3.0), "the bottom halo must come from the SW tile" def test_block_straddling_a_boundary_is_assembled_from_both_files(tmp_path): - """Um bloco que cai metade num tile e metade no outro precisa das duas - metades — é o caso que quebra se a montagem for feita tile a tile.""" + """A block that falls half in one tile and half in the other needs + both halves — the case that breaks if assembly is done tile by tile.""" ws = _ws(tmp_path, (2 * T, 2 * T), block=4, halo=1) - load_zarr_tiles_into_workspace(ws, _quadrantes(tmp_path, (1.0, 2.0, 3.0, 4.0))) - lido = MemmapRasterWorkspace(tmp_path / "ws") - nucleo = lido.read_block_core(Block(r0=0, r1=4, c0=T - 2, c1=T + 2), "v") - assert np.all(nucleo[:, :2] == 1.0) and np.all(nucleo[:, 2:] == 2.0) + load_zarr_tiles_into_workspace(ws, _quadrants(tmp_path, (1.0, 2.0, 3.0, 4.0))) + reopened = MemmapRasterWorkspace(tmp_path / "ws") + core_values = reopened.read_block_core(Block(r0=0, r1=4, c0=T - 2, c1=T + 2), "v") + assert np.all(core_values[:, :2] == 1.0) and np.all(core_values[:, 2:] == 2.0) -# ── buracos da malha ───────────────────────────────────────────────────────── +# ── holes in the tiling ────────────────────────────────────────────────────── def test_missing_tile_becomes_fill_not_garbage(tmp_path): - """Buraco real da malha (o MapBiomas tem vários) vira nodata.""" + """A real hole in the tiling (MapBiomas exports have several) becomes nodata.""" ws = _ws(tmp_path, (2 * T, 2 * T)) - load_zarr_tiles_into_workspace(ws, _quadrantes(tmp_path, (1.0, 2.0, None, 4.0))) - lido = MemmapRasterWorkspace(tmp_path / "ws") - buraco = lido.read_block_core(Block(r0=T, r1=T + 4, c0=0, c1=4), "v") - assert np.all(np.isnan(buraco)) + load_zarr_tiles_into_workspace(ws, _quadrants(tmp_path, (1.0, 2.0, None, 4.0))) + reopened = MemmapRasterWorkspace(tmp_path / "ws") + hole = reopened.read_block_core(Block(r0=T, r1=T + 4, c0=0, c1=4), "v") + assert np.all(np.isnan(hole)) def test_halo_over_a_hole_is_fill_while_core_keeps_data(tmp_path): ws = _ws(tmp_path, (2 * T, 2 * T), block=4, halo=1) - load_zarr_tiles_into_workspace(ws, _quadrantes(tmp_path, (1.0, 2.0, None, 4.0))) - lido = MemmapRasterWorkspace(tmp_path / "ws") - blk = Block(r0=T - 4, r1=T, c0=0, c1=4) # último bloco do tile NO - jan = lido.read_block_with_halo(blk, boundary_value=0.0)["v"] - assert np.all(jan[1:-1, 1:-1] == 1.0), "núcleo mantém o dado" - assert np.all(np.isnan(jan[-1, 1:-1])), "halo cai no buraco -> fill" + load_zarr_tiles_into_workspace(ws, _quadrants(tmp_path, (1.0, 2.0, None, 4.0))) + reopened = MemmapRasterWorkspace(tmp_path / "ws") + blk = Block(r0=T - 4, r1=T, c0=0, c1=4) # last block of the NW tile + win = reopened.read_block_with_halo(blk, boundary_value=0.0)["v"] + assert np.all(win[1:-1, 1:-1] == 1.0), "the core keeps the data" + assert np.all(np.isnan(win[-1, 1:-1])), "the halo falls in the hole -> fill" def test_explicit_fill_value_is_used(tmp_path): ws = _ws(tmp_path, (2 * T, 2 * T)) load_zarr_tiles_into_workspace( - ws, _quadrantes(tmp_path, (1.0, 2.0, None, 4.0)), fill=-9999.0 + ws, _quadrants(tmp_path, (1.0, 2.0, None, 4.0)), fill=-9999.0 ) - lido = MemmapRasterWorkspace(tmp_path / "ws") - assert np.all(lido.read_block_core(Block(r0=T, r1=T + 4, c0=0, c1=4), "v") == -9999.0) + reopened = MemmapRasterWorkspace(tmp_path / "ws") + assert np.all(reopened.read_block_core(Block(r0=T, r1=T + 4, c0=0, c1=4), "v") == -9999.0) def test_skip_empty_blocks_leaves_them_zero(tmp_path): - """Modo de economia de disco: bloco vazio não é escrito, então lê zero - (não `fill`) — a troca documentada no README.""" + """Disk-saving mode: an empty block is not written, so it reads zero + (not `fill`) — the trade-off documented in the README.""" ws = _ws(tmp_path, (2 * T, 2 * T)) load_zarr_tiles_into_workspace( - ws, _quadrantes(tmp_path, (1.0, 2.0, None, 4.0)), skip_empty_blocks=True + ws, _quadrants(tmp_path, (1.0, 2.0, None, 4.0)), skip_empty_blocks=True ) - lido = MemmapRasterWorkspace(tmp_path / "ws") - assert np.all(lido.read_block_core(Block(r0=T, r1=T + 4, c0=0, c1=4), "v") == 0.0) + reopened = MemmapRasterWorkspace(tmp_path / "ws") + assert np.all(reopened.read_block_core(Block(r0=T, r1=T + 4, c0=0, c1=4), "v") == 0.0) -# ── contrato e erros ───────────────────────────────────────────────────────── +# ── contract and errors ────────────────────────────────────────────────────── def test_single_tile_covering_the_grid_also_works(tmp_path): - """Variável global (sem tiles) chega como layout de um item só.""" + """A global variable (no tiles) arrives as a one-item layout.""" ws = _ws(tmp_path, (T, T)) url = _write_tile(tmp_path / "g.zarr", np.arange(T * T, dtype="float64").reshape(T, T)) load_zarr_tiles_into_workspace(ws, [{ "tile_id": None, "variable": "v", "url": url, "row_off": 0, "col_off": 0, "height": T, "width": T, }]) - lido = MemmapRasterWorkspace(tmp_path / "ws") - assert lido.read_block_core(Block(r0=0, r1=2, c0=0, c1=2), "v")[0, 0] == 0.0 + reopened = MemmapRasterWorkspace(tmp_path / "ws") + assert reopened.read_block_core(Block(r0=0, r1=2, c0=0, c1=2), "v")[0, 0] == 0.0 def test_array_name_can_differ_from_variable_name(tmp_path): @@ -201,7 +203,7 @@ def test_array_name_can_differ_from_variable_name(tmp_path): def test_empty_tile_list_raises(tmp_path): ws = _ws(tmp_path, (T, T)) - with pytest.raises(ValueError, match="vazia"): + with pytest.raises(ValueError, match="empty"): load_zarr_tiles_into_workspace(ws, []) @@ -214,30 +216,30 @@ def test_missing_key_names_what_is_missing(tmp_path): def test_undeclared_array_raises(tmp_path): ws = _ws(tmp_path, (T, T)) url = _write_tile(tmp_path / "g.zarr", np.ones((T, T))) - with pytest.raises(ValueError, match="não declara"): + with pytest.raises(ValueError, match="does not declare"): load_zarr_tiles_into_workspace(ws, [{ - "tile_id": None, "variable": "inexistente", "url": url, + "tile_id": None, "variable": "nonexistent", "url": url, "row_off": 0, "col_off": 0, "height": T, "width": T, }]) def test_tile_outside_workspace_raises(tmp_path): - """Um tile fora dos limites significa grade errada — falhar alto evita - um mosaico truncado que passaria despercebido.""" + """A tile out of bounds means a wrong grid — failing loudly avoids a + truncated mosaic that would go unnoticed.""" ws = _ws(tmp_path, (T, T)) url = _write_tile(tmp_path / "g.zarr", np.ones((T, T))) - with pytest.raises(ValueError, match="não cabe"): + with pytest.raises(ValueError, match="does not fit"): load_zarr_tiles_into_workspace(ws, [{ - "tile_id": "fora", "variable": "v", "url": url, + "tile_id": "outside", "variable": "v", "url": url, "row_off": T, "col_off": 0, "height": T, "width": T, }]) -# ── o contrato com o disscube, quando ele está disponível ──────────────────── +# ── the contract with disscube, when it is available ──────────────────────── def test_layout_shape_matches_disscube(tmp_path): - """Fixa que as chaves que este loader exige são as que o - CubeClient.tile_layout() produz. Sem o disscube instalado, pula.""" + """Pins that the keys this loader requires are the ones + CubeClient.tile_layout() produces. Skipped without disscube.""" pytest.importorskip("disscube") from disscube.client import CubeClient from disscube.models import DerivedVariable, GridSpec, SpatialSource @@ -256,23 +258,23 @@ def test_layout_shape_matches_disscube(tmp_path): derivation_id="d", spec_hash="h", tile_id="T1", asset_url="x.zarr", )) - exigidas = {"url", "variable", "row_off", "col_off", "height", "width"} - assert exigidas <= set(cube.tile_layout("v", "G")[0]) + required = {"url", "variable", "row_off", "col_off", "height", "width"} + assert required <= set(cube.tile_layout("v", "G")[0]) -# ── posições sobrepostas ───────────────────────────────────────────────────── -# Defesa em profundidade: mesmo que o layout venha errado de qualquer origem, -# dois pedaços na mesma posição não devem ser aceitos. O caso real que motivou -# isto: uma variável temporal cujo layout misturava as fatias de vários anos, -# todas nas mesmas posições — carregar isso deixaria o último ano vencer, sem -# erro nenhum. +# ── overlapping positions ──────────────────────────────────────────────────── +# Defence in depth: even if the layout comes out wrong from any source, two +# pieces at the same position must not be accepted. The real case behind +# this: a temporal variable whose layout mixed the slices of several years, +# all at the same positions — loading it would let the last year win, with +# no error at all. def test_two_tiles_at_the_same_position_raise(tmp_path): ws = _ws(tmp_path, (T, T)) a = _write_tile(tmp_path / "a.zarr", np.ones((T, T))) b = _write_tile(tmp_path / "b.zarr", np.full((T, T), 2.0)) base = {"variable": "v", "row_off": 0, "col_off": 0, "height": T, "width": T} - with pytest.raises(ValueError, match="ocupam a posição"): + with pytest.raises(ValueError, match="occupy position"): load_zarr_tiles_into_workspace(ws, [ {**base, "tile_id": "1985", "url": a}, {**base, "tile_id": "1995", "url": b}, @@ -285,13 +287,13 @@ def test_overlap_error_names_both_tiles(tmp_path): base = {"variable": "v", "row_off": 0, "col_off": 0, "height": T, "width": T} with pytest.raises(ValueError) as exc: load_zarr_tiles_into_workspace(ws, [ - {**base, "tile_id": "primeiro", "url": a}, - {**base, "tile_id": "segundo", "url": a}, + {**base, "tile_id": "first", "url": a}, + {**base, "tile_id": "second", "url": a}, ]) - assert "primeiro" in str(exc.value) and "segundo" in str(exc.value) + assert "first" in str(exc.value) and "second" in str(exc.value) def test_distinct_positions_still_accepted(tmp_path): - """Guarda: a checagem não pode recusar um mosaico legítimo.""" + """Guard: the check must not reject a legitimate mosaic.""" ws = _ws(tmp_path, (2 * T, 2 * T)) - load_zarr_tiles_into_workspace(ws, _quadrantes(tmp_path, (1.0, 2.0, 3.0, 4.0))) + load_zarr_tiles_into_workspace(ws, _quadrants(tmp_path, (1.0, 2.0, 3.0, 4.0)))