From 8f6b79a3d44e4d852c0495aeefa56ede3b7d5a97 Mon Sep 17 00:00:00 2001 From: Sergio Souza Costa Date: Wed, 23 Sep 2026 08:30:43 -0300 Subject: [PATCH] improved components code --- docs/api.md | 2 +- docs/decisions.md | 226 ++++++++++++++++++ pyproject.toml | 5 + src/disslucc/components/allocation/clue.py | 27 ++- src/disslucc/components/allocation/clue_s.py | 17 +- src/disslucc/components/potential/linear.py | 7 + src/disslucc/components/potential/logistic.py | 13 +- src/disslucc/schemas.py | 12 +- 8 files changed, 293 insertions(+), 16 deletions(-) diff --git a/docs/api.md b/docs/api.md index cc0e19d..8b454df 100644 --- a/docs/api.md +++ b/docs/api.md @@ -77,7 +77,7 @@ PotentialLinearRegression( potential_data: list[list[RegressionSpec]], # [region][class] demand, # object with DemandProtocol land_use_types: list[str], - land_use_no_data: str | None = None, # class to exclude from the calculation (e.g. water) + land_use_no_data: str | None = None, # array in [0,1]: scales potential down by (1 - value); NOT a binary exclusion mask region_attr: str = "region", # region array in the backend; created as 1 if absent ) ``` diff --git a/docs/decisions.md b/docs/decisions.md index ccd9c3e..a340730 100644 --- a/docs/decisions.md +++ b/docs/decisions.md @@ -755,3 +755,229 @@ received when TerraME/LuccME's own active development moved on. Not yet done as of this entry; this repo's own migration status section in `README.md` already anticipates it ("will be released, tagged, archived, and kept citable once this migration completes"). + +## Readability pass on `potential/`, for PIBIC onboarding: `scipy.special.expit` adopted, `numpy.ma` and `tensordot` rejected + +Requested while preparing onboarding material for PIBIC students who +will read `potential/linear.py` and `potential/logistic.py`: "does +anything here have an obvious simplification using something the +Python/NumPy/SciPy community already provides?" Three candidates were +tried; only one was kept, for reasons worth recording so they aren't +re-litigated later. + +**Adopted -- `logistic.py`'s hand-rolled sigmoid replaced with +`scipy.special.expit`.** The code manually clipped `z` to +`[-50, 50]` before `1 / (1 + exp(-z))`, with a comment explaining the +clip exists only to avoid `exp()` overflow on nodata cells. `expit` +is SciPy's numerically stable sigmoid -- exactly what that comment +was working around by hand -- so the clip and the manual formula both +go away. `scipy` was already an installed, transitive dependency (via +`dissmodel`'s own deps), but was not a *declared* one; added it +directly to `pyproject.toml`'s `dependencies`, since `disslucc` now +imports it itself rather than relying on an accident of what +`dissmodel` happens to pull in. Verified bit-for-bit identical output +on both Lab1 and Lab15 (same `mae`/`rmse`/`quantity_disagreement`/ +`allocation_disagreement` to full float precision, not just "within +tolerance") before and after -- `docs/validation.md` needs no update. + +**Rejected -- `numpy.ma` for `land_use_no_data` in `linear.py`.** The +suggestion was to replace `reg = reg * (1.0 - no_data_arr)` with a +masked array, on the reasoning that `land_use_no_data` sounds like a +binary "missing data" flag. Checked the actual values behind it first +(`data/input/csAC.zip`'s `outros` column, as loaded for the Lab1 +scenario): real-valued, 2183 distinct values between ~0 and 1, not a +0/1 mask. It's the *fraction* of a cell already committed to a +non-transitionable "other" class, and potential is scaled down +proportionally -- not zeroed out. `numpy.ma` is for hard/boolean +masking; using it here would have silently changed the science despite +looking like a pure refactor. Left the line as-is, added a comment +recording this so the next person doesn't repeat the same +name-implies-semantics mistake. Good example, for the PIBIC material +itself, of why to check the data before trusting a parameter's name. + +**Rejected (for now) -- vectorizing `const + sum(beta_k * driver_k)` +with `np.stack`/`np.tensordot` instead of the Python `for col, beta in +spec.betas.items()` loop.** Mathematically the same linear +combination, and arguably more idiomatic NumPy for a small number of +drivers. Tried it in both `linear.py` and `logistic.py` and reran the +full validation: Lab15 (binary, exact-match assertions) was +unaffected, but Lab1's `mae` moved from `0.0035832335619405574` to +`0.003583233421432411` -- floating-point summation reordering, not a +bug, and still far inside the `0.01` tolerance both +`tests/test_validation_lab1.py` and `disslucc-continuous`'s own +criterion use. Reverted anyway: this repo's own rule (see top of this +file, and `docs/validation.md`'s header) is that `docs/validation.md`'s +numbers are cited in `dissmodel`'s JOSS paper, and "the rounded number +published doesn't change" is a weaker bar than "bit-identical unless a +change is deliberately about the science." A cosmetic vectorization +isn't worth even a last-decimal-digit drift here. Revisit only if a +real performance need shows up (many drivers, large rasters) that +would justify accepting that tradeoff explicitly and re-publishing +`docs/validation.md`'s numbers alongside it. + +## Same pass, `allocation/clue.py`: two safe cosmetic fixes, one documented-not-fixed dead branch + +Continuation of the readability pass above, this time over +`AllocationClueLike` (the CLUE-like continuous allocator). Same +discipline: verify bit-for-bit identical Lab1/Lab15 output, run the +full suite, before keeping any change. + +**Applied -- two zero-risk stylistic fixes in `_correct_cell_change`.** +`flat = lambda lu: ...` (a lambda assigned to a name -- PEP 8 E731, and +it shows as `` in tracebacks instead of a real function name) +became a nested `def flat(lu: str) -> np.ndarray: ...`. Needed a +`cast(np.ndarray, ...)` inside it to keep mypy clean (`self.backend.get()` +returns `Any`), matching the pattern `_mask()` already uses just above +it. `[self.allocation_data[i].min_value for i in range(len(lus))]` (and +the `max_value` line next to it) became +`[spec.min_value for spec in self.allocation_data]` -- `allocation_data` +is already in `land_use_types` order, so indexing through `range(len(...))` +bought nothing. Verified identical Lab1 `mae`/`rmse` and Lab15 +disagreement to full float precision before/after; `pytest`/`mypy`/`ruff` +all clean. + +**Found, documented, deliberately not fixed -- the `no_data` exclusion +in `_apply_complementar`'s deficit-correction branch is dead on two +levels.** `no_data = getattr(self, "land_use_no_data", None)`: that +attribute is a `PotentialLinearRegression.setup()` parameter (see +`potential/linear.py`) that `AllocationClueLike.setup()` never accepts +-- so `no_data` is always `None`, and `lu != None` is always `True` for +a class name, meaning the `eligible = [lu for lu in others if lu != +no_data]` filter never actually excludes anything. Traced the likely +original intent by checking Lab1's real config +(`examples/run_lab1_real.py`): `land_use_no_data="outros"` (on +`Potential`) and `static["outros"] == 1` (on `Allocation`) name the +same class, so this was almost certainly meant to protect the +region's fixed/static class from this branch's deficit correction -- +the same protection `_compute_change`/`_correct_cell_change` already +give `static == 1` classes elsewhere in this same file, just spelled +with the wrong variable here. Also instrumented all 7 Lab1 steps: +`deficit.any()` was never `True` once -- this whole branch is +unexercised by every scenario in the repository, with no test proving +either the current or a fixed version behaves correctly. + +Decided not to fix this now, on both counts: (1) the correct fix +isn't "wire `land_use_no_data` through" -- it's changing the criterion +to `self.static[lu] != 1`, a different and more invasive change than +it looks like at first; (2) there is zero test coverage of this +branch, current or fixed, and this repository's whole discipline +(`docs/validation.md` numbers cited in `dissmodel`'s JOSS paper) is +built on not changing unverified behavior without a benchmark to +check it against. A synthetic scenario that forces `deficit.any()` +would be the right way to close this gap, but that's new test-writing +work, not a readability pass. Recorded here, and as an inline comment +at the call site, so the next person doesn't need to re-derive this +from scratch, and doesn't mistake "wiring `land_use_no_data` through" +for the actual fix. + +## Same pass, `allocation/clue_s.py`: one applied, two considered and rejected + +Continuation of the readability pass, now over `AllocationDClueSLike` +(discrete CLUE-S allocation). Same discipline as the two entries above: +bit-for-bit Lab1/Lab15 comparison before keeping anything. + +**Applied.** `max_diff = float(np.max(np.abs(list(diff.values()))))` -- +`diff` is a plain `dict[str, float]` with one scalar per land-use class +(a handful of classes, not per-cell). Round-tripping a handful of +Python floats through `list -> np.array -> np.abs -> np.max -> float` +buys nothing; `max(abs(v) for v in diff.values())` is the same value, +plain stdlib. Verified bit-for-bit identical Lab1/Lab15 metrics; +`pytest`/`mypy`/`ruff` all clean. + +**Considered, rejected -- `numpy.ma` for the masked argmax +(`scores = np.where(allowed, scores, -np.inf); np.argmax(scores, axis=1)`).** +This *looks* like the textbook `numpy.ma` use case (a real boolean mask, +unlike the `land_use_no_data` false alarm in `clue.py`), but the +`-np.inf` sentinel is already the standard, well-understood idiom for +"argmax over allowed choices only", and masked-array argmax has its own +sharp edge -- behavior on a fully-masked row (a cell with zero allowed +transitions) isn't obviously the same as `-np.inf`'s behavior (picks +index 0 deterministically) without checking. No test exercises a +transition matrix with a fully-blocked cell, so there's no way to +verify equivalence. Not worth the swap for a purely cosmetic gain over +already-correct, already-idiomatic code. + +**Considered, rejected -- restructuring the `for n_iter in +range(self.max_iteration + 1): ... if n_iter >= self.max_iteration: +raise` convergence loop** into a cleaner `for/else`. Mechanically +sound, but touches the exact iteration count at which non-convergence +raises `RuntimeError`, and no test in this repository exercises that +failure path (Lab1/Lab15 always converge well within their configured +`max_iteration`). Changing loop bounds with zero test coverage of the +boundary condition it changes is exactly the kind of "looks safe, +unverifiable" edit this readability pass has been avoiding throughout +`potential/` and `allocation/`. + +## `allocation/clue_s.py`: hoisting `pot`/`1+tau` out of the convergence loop -- an actual algorithmic improvement, not just readability + +Different category from the three entries above: this one changes +*what the algorithm does* (fewer redundant recomputations per call), +not just how it's written, so it got measured, not just diffed. + +`pot = np.stack([self.backend.get(lu + "_pot").ravel() for lu in +lu_types], axis=1)` was being recomputed **every iteration** of +`execute()`'s convergence loop, even though `_pot` is written once +by `Potential.execute()` before `Allocation.execute()` runs at all, and +nothing inside this loop ever writes it back -- only `iter_vec` changes +between iterations. Same for `1.0 + tau`. Instrumented Lab15 first +(before touching anything): `execute()` took between 1 and 68 +iterations per time step across the 6 steps, so this was a real, +measurable amount of repeated work, not a one-off. + +Hoisted both `pot` and `one_plus_tau = 1.0 + tau` above the loop. +Verified bit-for-bit identical Lab1/Lab15 metrics (this change only +touches `clue_s.py`, so Lab1 was never going to move, but checked +anyway per this session's own discipline). Measured wall-clock time +for `run_lab15_raster()` (8 repetitions each, sorted): + + before: median 0.4792s, min 0.4510s + after: median 0.4283s, min 0.4077s (~10-11% faster end to end) + +That end-to-end number includes shapefile loading (a fixed cost this +change doesn't touch), so the saving *inside* the allocation loop +itself is proportionally larger and should grow with raster size and +iteration count -- the avoided work scales with +`n_iterations * n_classes * n_cells` per step, while the I/O floor +stays constant. Not benchmarked at larger scale; if raster size grows +significantly in a future scenario, worth re-measuring rather than +assuming the same ~10% holds. + +`pytest`/`mypy`/`ruff` all clean after the change. + +## Ecosystem-wide review pass (PIBIC-focused): `RegressionSpec.newconst` footgun and an `api.md` inaccuracy + +Broader pass across the whole package (not just `potential/`/`allocation/` +this time), specifically looking for anything that would confuse an +undergraduate reading this code for the first time. Two real findings, +both fixed. + +**`schemas.RegressionSpec.newconst` was a silent-no-op constructor +trap.** It's a regular dataclass field (`newconst: float = 0.0`), +so `RegressionSpec(const=0.5, newconst=99)` was accepted -- but +`PotentialLinearRegression.setup()`/`execute()` always overwrite it +with `spec.newconst = spec.const` before ever reading it, on every +single step, not just the first. So any value passed at construction +time was silently discarded, immediately. Confirmed with `grep` that +nothing in `src/`/`examples/`/`tests/` ever constructs a +`RegressionSpec` with `newconst=...` -- this was a purely latent trap +for a future caller, not a bug that ever fired. Fixed with +`field(default=0.0, init=False, repr=False, compare=False)`: +`newconst` is internal runtime state the model manages, not something +a caller configures, so it's no longer part of `__init__`, doesn't +clutter `repr()`, and doesn't affect `==` between two specs that only +differ in runtime-mutated state. Verified: `RegressionSpec(const=0.5, +newconst=99)` now raises `TypeError` immediately instead of silently +accepting and discarding the value. Bit-for-bit identical Lab1/Lab15 +metrics; `pytest`/`mypy`/`ruff` all clean. + +**`docs/api.md` mischaracterized `land_use_no_data`.** It described the +parameter as "class to exclude from the calculation (e.g. water)" -- +which is exactly the wrong mental model this session already +disproved by checking the real data (see the `potential/linear.py` +entry above): it's a real-valued `[0,1]` array that scales potential +down proportionally (`reg * (1 - value)`), not a binary +exclusion/no-data mask. This is the same misconception that produced +the dead `no_data` code in `allocation/clue.py`, now also fixed in the +one document a student would read *before* the source, making it more +likely to mislead, not less. Corrected the parameter comment in +`api.md` to describe the actual behavior. diff --git a/pyproject.toml b/pyproject.toml index 999e8a8..98c084f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,6 +8,11 @@ license = { text = "MIT" } requires-python = ">=3.11" dependencies = [ "dissmodel>=0.6.4,<0.7.0", + # scipy.special.expit (potential/logistic.py) -- was already pulled + # in transitively by dissmodel's own deps, but disslucc now imports + # it directly, so it needs to be a declared dependency, not an + # accident of what dissmodel happens to require today. + "scipy", ] [project.optional-dependencies] diff --git a/src/disslucc/components/allocation/clue.py b/src/disslucc/components/allocation/clue.py index 96daec2..ca7b772 100644 --- a/src/disslucc/components/allocation/clue.py +++ b/src/disslucc/components/allocation/clue.py @@ -146,15 +146,17 @@ def _correct_cell_change(self) -> None: TOL = 0.005 MAX_L = 25 - alloc_min = np.array([self.allocation_data[i].min_value for i in range(len(lus))]) - alloc_max = np.array([self.allocation_data[i].max_value for i in range(len(lus))]) + alloc_min = np.array([spec.min_value for spec in self.allocation_data]) + alloc_max = np.array([spec.max_value for spec in self.allocation_data]) static = np.array([self.static[lu] for lu in lus]) mask = self._mask() shape = self.shape flat_mask = mask.ravel() - flat = lambda lu: self.backend.get(lu).ravel().astype(np.float32) + def flat(lu: str) -> np.ndarray: + return cast(np.ndarray, self.backend.get(lu).ravel().astype(np.float32)) + originals = {lu: flat(lu) for lu in lus} vals = np.stack([flat(lu) for lu in lus], axis=1) pasts = np.stack([flat(lu + "_past") for lu in lus], axis=1) @@ -238,6 +240,25 @@ def _apply_complementar(self) -> None: comp_original = self.backend.get(self.complementar_lu).astype(np.float32) self.backend.arrays[self.complementar_lu] = np.where(mask, comp, comp_original) + # NOTE: this whole branch is unexercised by every scenario in this + # repository -- instrumented Lab1's 7 steps and `deficit.any()` was + # never True. No test proves it behaves correctly, as-is or fixed. + # + # `no_data` is also dead on arrival: `land_use_no_data` is a + # PotentialLinearRegression.setup() parameter (see potential/linear.py) + # that AllocationClueLike.setup() never accepts, so + # `getattr(self, "land_use_no_data", None)` is always None here, and + # `lu != None` is always True -- the filter below never actually + # excludes anything. In Lab1, `land_use_no_data="outros"` and + # `static["outros"] == 1` happen to name the same class, which is + # almost certainly the real intent: protect the region's fixed/static + # class from being reduced just to close this deficit, the same way + # `_compute_change`/`_correct_cell_change` already protect it via + # `self.static[lu] == 1` / `is_static`. If this branch is ever made + # reachable (a synthetic scenario that forces `deficit.any()`), fix + # the criterion to `self.static[lu] != 1`, not to wiring + # `land_use_no_data` through -- that would still be the wrong + # protection, just no longer silently broken. deficit = np.maximum(0.0, total - 1.0) if deficit.any(): no_data = getattr(self, "land_use_no_data", None) diff --git a/src/disslucc/components/allocation/clue_s.py b/src/disslucc/components/allocation/clue_s.py index 122493e..03dc785 100644 --- a/src/disslucc/components/allocation/clue_s.py +++ b/src/disslucc/components/allocation/clue_s.py @@ -95,9 +95,17 @@ def execute(self) -> None: iter_vec = np.zeros(n_lu, dtype=np.float64) + # `_pot` is written once by Potential.execute(), before this + # component runs, and nothing in this loop ever writes it back -- + # only iter_vec changes between iterations. Hoisting both `pot` and + # `1.0 + tau` out of the loop skips re-reading/re-stacking backend + # arrays on every one of the (often dozens of) convergence + # iterations, for the same result every time. + pot = np.stack([self.backend.get(lu + "_pot").ravel() for lu in lu_types], axis=1) + one_plus_tau = 1.0 + tau + for n_iter in range(self.max_iteration + 1): - pot = np.stack([self.backend.get(lu + "_pot").ravel() for lu in lu_types], axis=1) - scores = (1.0 + tau) * pot + iter_vec[np.newaxis, :] + scores = one_plus_tau * pot + iter_vec[np.newaxis, :] scores = np.where(allowed, scores, -np.inf) best_lu_idx = np.argmax(scores, axis=1) # (n_cells,) @@ -107,7 +115,10 @@ def execute(self) -> None: self.backend.arrays[lu] = np.where(mask, new_arr, self.backend.get(lu)) diff = self._calc_diff(mask) - max_diff = float(np.max(np.abs(list(diff.values())))) + # diff has one scalar per land-use class (a handful) -- plain + # max()/abs() is the idiomatic way to reduce that, no need to + # round-trip through a NumPy array for a handful of Python floats. + max_diff = max(abs(v) for v in diff.values()) if max_diff <= self.max_difference: break diff --git a/src/disslucc/components/potential/linear.py b/src/disslucc/components/potential/linear.py index e33302b..49ea3c6 100644 --- a/src/disslucc/components/potential/linear.py +++ b/src/disslucc/components/potential/linear.py @@ -90,6 +90,13 @@ def _compute_potential(self, r_number: int, lu_idx: int) -> None: reg = np.clip(reg, 0.0, 1.0) if self.land_use_no_data: + # NOTE: despite the parameter's name, this is not a binary + # "missing data" mask -- it's the fraction of the cell already + # committed to a non-transitionable "other" class (real-valued, + # e.g. 0.37), so potential is scaled down proportionally, not + # zeroed out. Confirmed against data/input/csAC.zip's `outros` + # column before touching this line -- don't replace it with + # numpy.ma or any hard/boolean mask. no_data_arr = self.backend.get(self.land_use_no_data).astype(np.float32) reg = reg * (1.0 - no_data_arr) diff --git a/src/disslucc/components/potential/logistic.py b/src/disslucc/components/potential/logistic.py index af460e8..0b8e888 100644 --- a/src/disslucc/components/potential/logistic.py +++ b/src/disslucc/components/potential/logistic.py @@ -22,6 +22,7 @@ import numpy as np from dissmodel.geo import SyncRasterModel +from scipy.special import expit from ...schemas import LogisticRegressionSpec @@ -73,14 +74,10 @@ def _compute_potential(self, mask: np.ndarray, lu_idx: int, spec: LogisticRegres for col, beta in spec.betas.items(): z = z + beta * self.backend.get(col).astype(np.float32) - # cells outside the mask (filled with nodata, e.g. -1) can produce - # extreme z with the real coefficients -- clip avoids overflow in - # exp() without changing the result (sigmoid already saturates - # well before +-50) - z = np.clip(z, -50.0, 50.0) - - # numerically stable sigmoid - prob = 1.0 / (1.0 + np.exp(-z)) + # scipy.special.expit is the standard, numerically stable sigmoid + # (no manual overflow clipping needed -- it handles extreme z from + # cells outside the mask, e.g. nodata=-1, without over/underflow). + prob = expit(z) elas = np.where(self.backend.get(lu) == 1, spec.elasticity, 0.0).astype(np.float32) pot = prob + elas diff --git a/src/disslucc/schemas.py b/src/disslucc/schemas.py index ab556db..b0a5e26 100644 --- a/src/disslucc/schemas.py +++ b/src/disslucc/schemas.py @@ -21,11 +21,21 @@ class RegressionSpec: explicitly: it's the `const` value already adjusted by `_adapt_constants()` each step -- mutable by design, not a typo for `const`. + + `init=False`: `PotentialLinearRegression.setup()`/`execute()` + always overwrite it with `spec.newconst = spec.const` before ever + reading it (every step, not just once), so a value passed to the + constructor would be silently discarded on the first run -- this + is internal runtime state the model manages, not something a + caller configures. Excluded from `__init__`/`repr`/`==` so + `RegressionSpec(const=0.5, newconst=99)` is a clear TypeError + instead of a silent no-op, and two specs that only differ in + runtime-mutated state still compare equal. """ const: float betas: dict[str, float] = field(default_factory=dict) is_log: bool = False - newconst: float = 0.0 + newconst: float = field(default=0.0, init=False, repr=False, compare=False) @dataclass