diff --git a/CLAUDE.md b/CLAUDE.md index 12951dd..55469e2 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -59,7 +59,7 @@ The reference results live in `benchmark/goldens/`, a copy of the goldens generated in [LambdaGeo/terrame-docker](https://github.com/LambdaGeo/terrame-docker) v0.1.1, which keeps the generating scripts, the original TerraME outputs, the generator and goldens for all 21 LuccME labs. This repository keeps only the -goldens its tests use (`lab01`, `lab15`: the LuccME package's labs; +goldens its tests use (`lab01`, `lab15`, `lab03`, `lab06`: the LuccME package's labs; `lab01_md1643`, `lab15_md10`); add one together with the component and test that need it, never ahead of time. `lab01_md1643` and `lab15_md10` are the scenarios behind `docs/validation.md` (same coefficients and demand as `lab01`/`lab15`, but `maxDifference` 1643 @@ -81,7 +81,7 @@ PR/commit and say so explicitly -- don't let it drift silently. ```bash python -m venv venv && source venv/bin/activate pip install -e ".[examples,dev]" -pytest tests/ -v # expect: 19 passed, 2 xfailed +pytest tests/ -v # expect: 73 passed, 2 xfailed (38 passed, 1 skipped without lupa) mypy src/disslucc # expect: clean ``` diff --git a/benchmark/README.md b/benchmark/README.md index d42a2e2..7956ac5 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -30,6 +30,8 @@ implemented, together with the test that uses them. | `lab01_md1643` | same, `maxDifference` 1643: iterates up to 26 times per year | `test_goldens_per_year.py`, `test_validation_lab1.py`, discriminance | | `lab15` | PreComputedValues + DLogisticRegression + DClueSLike (`maxDifference` 300) | `test_goldens_per_year.py` | | `lab15_md10` | same, `maxDifference` 10: iterates 56–67 times per year | `test_goldens_per_year.py`, `test_validation_lab15.py`, discriminance | +| `lab03` | PreComputedValues + CSpatialLagRegression + CClueLikeSaturation (`maxDifference` 1643) | `test_spatial_lag_golden.py`, `test_saturation_golden.py` | +| `lab06` | same as `lab03`, plus `updateYears = {2009}` (`ti` from `csAC_2009`, `data/input/csAC_2009.zip`) | `test_spatial_lag_golden.py`, `test_saturation_golden.py` | The last year of `lab01_md1643` and of `lab15_md10` is this repository's former reference (`benchmark/data/*.zip`). diff --git a/benchmark/goldens/lab03/lab03.csv.gz b/benchmark/goldens/lab03/lab03.csv.gz new file mode 100644 index 0000000..b316365 Binary files /dev/null and b/benchmark/goldens/lab03/lab03.csv.gz differ diff --git a/benchmark/goldens/lab03/manifest.json b/benchmark/goldens/lab03/manifest.json new file mode 100644 index 0000000..2946a8d --- /dev/null +++ b/benchmark/goldens/lab03/manifest.json @@ -0,0 +1,45 @@ +{ + "lab": "lab03", + "source": { + "script": "luccme/tests/functional/lab03.lua", + "sha256": "f32e13c853c05b0ea3df9368086e3ccc4fab8d3ba8322f39f15a7f056160a027" + }, + "engine": { + "image": "terrame-luccme", + "terrame": "2.0.1", + "terralib": "5.5.1", + "luccme_commit": "6244dd461f94259efb6e1d2170d32fc7e28c033c" + }, + "generator": { + "script": "benchmark/harness.lua", + "sha256": "5f7504f11442d83981565676635ffd996034f721a5ad08ec0cab082d2e948a97" + }, + "status": "ok", + "years": [ + 2008, + 2014 + ], + "n_cells": 6574, + "columns": [ + "f_out", + "d_out", + "outros_out", + "f_pot", + "d_pot", + "outros_pot" + ], + "file": { + "name": "lab03.csv.gz", + "rows": 46018, + "sha256": "dba249299d94626100b97af14cc46734a17df2a957ed36bda9ff3881c6ebba29" + }, + "crosscheck_vs_original_output": { + "tolerance": 1e-09, + "max_abs_diff": { + "Lab03_2014.dbf:d_out": 5.000002148425331e-13 + }, + "note": "last year of the CSV compared with the .dbf files of the output folder: the script's own output and, for references, the original TerraME output" + }, + "iterations_per_year": {}, + "iterations_note": "continuous: LuccME's 'Number of iterations'; discrete: largest n in 'Iteration -> n' (0 = first pass accepted)" +} diff --git a/benchmark/goldens/lab03/terrame.log b/benchmark/goldens/lab03/terrame.log new file mode 100644 index 0000000..51ee38e --- /dev/null +++ b/benchmark/goldens/lab03/terrame.log @@ -0,0 +1,118 @@ + +Verifying Model parameters +Verifying Demand parameters +Verifying Potential parameters +Verifying Allocation parameters + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2008 Step: 0 +f area: 137878 Difference: 0 +d area: 19982 Difference: 0 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.10000000002027 dir: -1 const: 0.05266679 -> 0.05266679 0 1 +d elas: 0.10000000001401 dir: -1 const: 0.01431553 -> 0.01431553 0 1 +outros elas: 0.10000016027163 dir: -1 const: 0 -> 0 0 1 + +Demand allocated correctly in 2008 Number of iterations: 0 Maximum error: 0.01040035000733 +[harness] lab03: Lab03, years 2008..2014, columns f_out,d_out,outros_out,f_pot,d_pot,outros_pot + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2009 Step: 0 +f area: 137418 Difference: -204 +d area: 20334 Difference: 96 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.099851813713511 dir: -1 const: 0.052648226568119 -> 0.052648226568119 -1.8563431881255e-05 0.99998143656812 +d elas: 0.0995264862247 dir: 1 const: 0.014443615865131 -> 0.014443615865131 0.00012808586513093 1.0001280858651 +outros elas: 0.10000016027163 dir: 0 const: 0.0 -> 0.0 0.0 1.0 + +Demand allocated correctly in 2009 Number of iterations: 0 Maximum error: 203.9372570531 + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2010 Step: 0 +f area: 137079 Difference: -287 +d area: 20673 Difference: 179 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.099791090223736 dir: -1 const: 0.052629628612048 -> 0.052629628612048 -1.8597956070319e-05 0.99998140204393 +d elas: 0.099132190192513 dir: 1 const: 0.014570081884318 -> 0.014570081884318 0.00012646601918755 1.0001264660192 +outros elas: 0.10000016027163 dir: 0 const: 0.0 -> 0.0 0.0 1.0 + +Demand allocated correctly in 2010 Number of iterations: 0 Maximum error: 286.971568782 + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2011 Step: 0 +f area: 136752 Difference: -358 +d area: 21000 Difference: 250 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.099739024534728 dir: -1 const: 0.052610995995854 -> 0.052610995995854 -1.8632616194334e-05 0.99998136738381 +d elas: 0.0988079001206 dir: 1 const: 0.014694968507161 -> 0.014694968507161 0.00012488662284242 1.0001248866228 +outros elas: 0.10000016027163 dir: 0 const: 0.0 -> 0.0 0.0 1.0 + +Demand allocated correctly in 2011 Number of iterations: 0 Maximum error: 357.82429921019 + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2012 Step: 0 +f area: 136436 Difference: -388 +d area: 21316 Difference: 280 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.099716648902286 dir: -1 const: 0.052590163421247 -> 0.052590163421247 -2.0832574607325e-05 0.99997916742539 +d elas: 0.098684256490536 dir: 1 const: 0.014832621312254 -> 0.014832621312254 0.00013765280509327 1.0001376528051 +outros elas: 0.10000016027163 dir: 0 const: 0.0 -> 0.0 0.0 1.0 + +Demand allocated correctly in 2012 Number of iterations: 0 Maximum error: 387.69424774163 + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2013 Step: 0 +f area: 136131 Difference: -408 +d area: 21622 Difference: 300 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.099701765439119 dir: -1 const: 0.052569287356422 -> 0.052569287356422 -2.0876064825124e-05 0.99997912393517 +d elas: 0.098610855076548 dir: 1 const: 0.014968405021323 -> 0.014968405021323 0.00013578370906851 1.0001357837091 +outros elas: 0.10000016027163 dir: 0 const: 0.0 -> 0.0 0.0 1.0 + +Demand allocated correctly in 2013 Number of iterations: 0 Maximum error: 407.20663381234 + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2014 Step: 0 +f area: 135836 Difference: -418 +d area: 21918 Difference: 310 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.09969381618728 dir: -1 const: 0.052548367612094 -> 0.052548367612094 -2.0919744327618e-05 0.99997908025567 +d elas: 0.098582216155697 dir: 1 const: 0.015102369708185 -> 0.015102369708185 0.00013396468686188 1.0001339646869 +outros elas: 0.10000016027163 dir: 0 const: 0.0 -> 0.0 0.0 1.0 + +Demand allocated correctly in 2014 Number of iterations: 0 Maximum error: 417.18589488399 + +Saving Lab03_2014. +Elapsed time: 00:00:04 hh:mm:ss + +End of Simulation +[harness] lab03: 5 files of the lab's own output in /work/out/lab03 +[harness] lab03: CSV written to /work/out/lab03 diff --git a/benchmark/goldens/lab06/lab06.csv.gz b/benchmark/goldens/lab06/lab06.csv.gz new file mode 100644 index 0000000..4175177 Binary files /dev/null and b/benchmark/goldens/lab06/lab06.csv.gz differ diff --git a/benchmark/goldens/lab06/manifest.json b/benchmark/goldens/lab06/manifest.json new file mode 100644 index 0000000..a48e5e6 --- /dev/null +++ b/benchmark/goldens/lab06/manifest.json @@ -0,0 +1,45 @@ +{ + "lab": "lab06", + "source": { + "script": "luccme/tests/functional/lab06.lua", + "sha256": "1248ddfea1df91ef6bd07aea3ff325fc748e10107883c548d69ae35a20041b5b" + }, + "engine": { + "image": "terrame-luccme", + "terrame": "2.0.1", + "terralib": "5.5.1", + "luccme_commit": "6244dd461f94259efb6e1d2170d32fc7e28c033c" + }, + "generator": { + "script": "benchmark/harness.lua", + "sha256": "5f7504f11442d83981565676635ffd996034f721a5ad08ec0cab082d2e948a97" + }, + "status": "ok", + "years": [ + 2008, + 2014 + ], + "n_cells": 6574, + "columns": [ + "f_out", + "d_out", + "outros_out", + "f_pot", + "d_pot", + "outros_pot" + ], + "file": { + "name": "lab06.csv.gz", + "rows": 46018, + "sha256": "456b5346f631028b7dae3e2cad247760e59abc57de40d549e02f2331f163d1b6" + }, + "crosscheck_vs_original_output": { + "tolerance": 1e-09, + "max_abs_diff": { + "Lab06_2014.dbf:d_out": 5.000002148425331e-13 + }, + "note": "last year of the CSV compared with the .dbf files of the output folder: the script's own output and, for references, the original TerraME output" + }, + "iterations_per_year": {}, + "iterations_note": "continuous: LuccME's 'Number of iterations'; discrete: largest n in 'Iteration -> n' (0 = first pass accepted)" +} diff --git a/benchmark/goldens/lab06/terrame.log b/benchmark/goldens/lab06/terrame.log new file mode 100644 index 0000000..9c55404 --- /dev/null +++ b/benchmark/goldens/lab06/terrame.log @@ -0,0 +1,121 @@ + +Verifying Model parameters +Verifying Demand parameters +Verifying Potential parameters +Verifying Allocation parameters + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2008 Step: 0 +f area: 137878 Difference: 0 +d area: 19982 Difference: 0 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.10000000002027 dir: -1 const: 0.05266679 -> 0.05266679 0 1 +d elas: 0.10000000001401 dir: -1 const: 0.01431553 -> 0.01431553 0 1 +outros elas: 0.10000016027163 dir: -1 const: 0 -> 0 0 1 + +Demand allocated correctly in 2008 Number of iterations: 0 Maximum error: 0.01040035000733 +[harness] lab06: Lab06, years 2008..2014, columns f_out,d_out,outros_out,f_pot,d_pot,outros_pot + +Updating dynamic variables... +layer_2009 + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2009 Step: 0 +f area: 137389 Difference: -233 +d area: 20331 Difference: 93 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.09983074185143 dir: -1 const: 0.052648226568119 -> 0.052648226568119 -1.8563431881255e-05 0.99998143656812 +d elas: 0.099542216292271 dir: 1 const: 0.014443615865131 -> 0.014443615865131 0.00012808586513093 1.0001280858651 +outros elas: 0.10000016027163 dir: 0 const: 0.0 -> 0.0 0.0 1.0 + +Demand allocated correctly in 2009 Number of iterations: 0 Maximum error: 232.93682142338 + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2010 Step: 0 +f area: 137053 Difference: -313 +d area: 20667 Difference: 173 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.099772292812196 dir: -1 const: 0.052629628612048 -> 0.052629628612048 -1.8597956070319e-05 0.99998140204393 +d elas: 0.099162815645913 dir: 1 const: 0.014570081884318 -> 0.014570081884318 0.00012646601918755 1.0001264660192 +outros elas: 0.10000016027163 dir: 0 const: 0.0 -> 0.0 0.0 1.0 + +Demand allocated correctly in 2010 Number of iterations: 0 Maximum error: 312.79287200209 + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2011 Step: 0 +f area: 136729 Difference: -381 +d area: 20991 Difference: 240 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.099722588869631 dir: -1 const: 0.052610995995854 -> 0.052610995995854 -1.8632616194334e-05 0.99998136738381 +d elas: 0.098851979284461 dir: 1 const: 0.014694968507161 -> 0.014694968507161 0.00012488662284242 1.0001248866228 +outros elas: 0.10000016027163 dir: 0 const: 0.0 -> 0.0 0.0 1.0 + +Demand allocated correctly in 2011 Number of iterations: 0 Maximum error: 380.35929244847 + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2012 Step: 0 +f area: 136417 Difference: -408 +d area: 21304 Difference: 268 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.099702461272719 dir: -1 const: 0.052590163421247 -> 0.052590163421247 -2.0832574607325e-05 0.99997916742539 +d elas: 0.098740566525434 dir: 1 const: 0.014832621312254 -> 0.014832621312254 0.00013765280509327 1.0001376528051 +outros elas: 0.10000016027163 dir: 0 const: 0.0 -> 0.0 0.0 1.0 + +Demand allocated correctly in 2012 Number of iterations: 0 Maximum error: 407.10642724726 + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2013 Step: 0 +f area: 136115 Difference: -424 +d area: 21607 Difference: 285 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.099689657048882 dir: -1 const: 0.052569287356422 -> 0.052569287356422 -2.0876064825124e-05 0.99997912393517 +d elas: 0.098678829724857 dir: 1 const: 0.014968405021323 -> 0.014968405021323 0.00013578370906851 1.0001357837091 +outros elas: 0.10000016027163 dir: 0 const: 0.0 -> 0.0 0.0 1.0 + +Demand allocated correctly in 2013 Number of iterations: 0 Maximum error: 423.73931471611 + +Executing Demand component +Executing Potential component +Executing Allocation component + +Year: 2014 Step: 0 +f area: 135822 Difference: -431 +d area: 21900 Difference: 293 +outros area: 6489 Difference: 0 +Region 1 +f elas: 0.099683763142858 dir: -1 const: 0.052548367612094 -> 0.052548367612094 -2.0919744327618e-05 0.99997908025567 +d elas: 0.098660722562989 dir: 1 const: 0.015102369708185 -> 0.015102369708185 0.00013396468686188 1.0001339646869 +outros elas: 0.10000016027163 dir: 0 const: 0.0 -> 0.0 0.0 1.0 + +Demand allocated correctly in 2014 Number of iterations: 0 Maximum error: 430.88351101935 + +Saving Lab06_2014. +Elapsed time: 00:00:05 hh:mm:ss + +End of Simulation +[harness] lab06: 5 files of the lab's own output in /work/out/lab06 +[harness] lab06: CSV written to /work/out/lab06 diff --git a/data/input/csAC_2009.zip b/data/input/csAC_2009.zip new file mode 100644 index 0000000..93634fc Binary files /dev/null and b/data/input/csAC_2009.zip differ diff --git a/docs/api.md b/docs/api.md index d2e9b06..eee8d44 100644 --- a/docs/api.md +++ b/docs/api.md @@ -117,6 +117,34 @@ feedback, it adjusts a global correction vector in --- +### `PotentialSpatialLagRegression` (continuous) + +Port of LuccME's `PotentialCSpatialLagRegression`, the potential of LuccME-BR. +For each class: `reg = newconst + Σ beta·x + ro·Y`, where `Y` is the mean share +of the class over the cell and its Moore neighbours (each share divided by +`1 - no-data share`; all-no-data neighbours left out); cells whose +`const + Σ beta·x + ro·Y` exceeds `max_reg` get 1, below `min_reg` get 0; then +`reg · (1 - no-data share)`, and `pot = reg - past share`. Writes `_pot` and +`_reg`. + +```python +PotentialSpatialLagRegression( + backend, + potential_data: list[list[SpatialLagRegressionSpec]], # [region][class] + demand, # DemandProtocol + land_use_types: list[str], + land_use_no_data: str | None = None, # share of the no-data class, as in PotentialLinearRegression + region_attr: str = "region", # created as 1 if absent + mask_attr: str = "mask", # cells that exist; created as 1 if absent +) +``` + +Unlike `PotentialLinearRegression`, the constant adaptation **accumulates** +(LuccME writes the adapted value back into `const` every step), and +`modify_driver(attr, rate)` ports LuccME's `modifyDriver` (used by +`AllocationClueLikeSaturation`). The pure function `spatial_lag_regression()` +in the same module computes one class, for tests and reuse. + ## Allocation ### `AllocationClueLike` (continuous) @@ -194,6 +222,47 @@ accepted). It matches TerraME in every year of `lab15_md10` (0, 67, 56, 56, 61, --- +### `AllocationClueLikeSaturation` (continuous) + +Port of LuccME's `AllocationCClueLikeSaturation`, the allocation of LuccME-BR: +`AllocationClueLike`'s CLUE loop plus a per-cell **saturation indicator** — +the share of the available area (not no-data, not protected) no longer in +`complementar_lu`, averaged over the 3 × 3 window with the cell, recomputed each +step. Where it exceeds a class's `change_limiar_value`, the change in the +demand's direction is halved or capped at `max_change_above_limiar`. + +```python +AllocationClueLikeSaturation( + backend, + demand, # DemandProtocol + potential, # RegionalPotentialProtocol + land_use_types: list[str], + allocation_data: list[list[SaturationAllocationSpec]], # [region][class], static per class + complementar_lu: str, + cell_area: float, + land_use_no_data: str | None = None, + attr_protection: str | None = None, # protected share, left out of the indicator + saturation_indicator: str = "saturationLimiar", # array written each step + max_difference: float = 1643, + max_iteration: int = 1000, + initial_elasticity: float = 0.1, + min_elasticity: float = 0.001, + max_elasticity: float = 1.5, + region_attr: str = "regionAloc", + mask_attr: str = "mask", + order_attr: str | None = None, # LuccME's cell order (see below); default row-major +) +``` + +> Unlike `AllocationClueLike`, this is LuccME's own `correctCellChange`, which +> runs in the Saturation variant (the region test is spelled right there). Its +> `BACKP` is declared outside the loop over cells, so a cell that lowers it lowers +> it for the cells visited after it: the result depends on the order of the +> cells, which `order_attr` gives (TerraME visits them in the order of the layer). +> +> `iterations_per_step` and `max_error_per_step` record, per step, what LuccME +> logs as "Number of iterations" and "Maximum error". + ## Validation ### `pontius_millones(pred, ref) -> dict` @@ -256,6 +325,25 @@ class AllocationSpec: max_value: float = 1.0 min_change: float = 0.0 max_change: float = 1.0 + +@dataclass +class SpatialLagRegressionSpec: + const: float # adapted every step, cumulatively + ro: float # spatial autoregressive coefficient + betas: dict[str, float] = {} + is_log: bool = False + min_reg: float = 0.0 + max_reg: float = 1.0 + +@dataclass +class SaturationAllocationSpec: + static: int = -1 + min_value: float = 0.0 + max_value: float = 1.0 + min_change: float = 0.0 + max_change: float = 1.0 + change_limiar_value: float = 1.0 # saturation above this limits the change... + max_change_above_limiar: float = 0.0 # ...to half, or to this ``` --- @@ -274,6 +362,10 @@ class DemandProtocol(Protocol): class PotentialProtocol(Protocol): def modify(self, r_number: int, lu_idx: int, direction: int) -> None: ... + +class RegionalPotentialProtocol(PotentialProtocol, Protocol): # AllocationClueLikeSaturation + potential_data: list # one entry per region + def modify_driver(self, attr_protection: str, rate: float) -> None: ... ``` --- @@ -299,9 +391,13 @@ premature abstraction. Provides: - `validate()` -- checks `record.source.uri` and `self.required_parameters` (list of required keys, declared per subclass) -- `load()` -- loads the GeoDataFrame via `load_dataset`, applies - `column_map`, **and already rasterizes** (`vector_to_raster_backend`) - -- returns the `RasterBackend` ready to use, not the GeoDataFrame. +- `load()` -- a `.tif`/`.tiff` source is a **GeoTIFF with named bands** + (a `name` tag per band, as `dissmodel.io.save_geotiff` writes): read as + it is, with its georeference, after checking that every land use and + driver has a band. Any other source is loaded as a GeoDataFrame via + `load_dataset`, gets `column_map`, **and is already rasterized** + (`vector_to_raster_backend`). Either way it returns the `RasterBackend` + ready to use, not the GeoDataFrame. Same convention as the real `LUCCRasterExecutor`, checked against the source: rasterizing is expensive, it runs once inside `load()`, never inside `run()`. @@ -310,8 +406,8 @@ premature abstraction. Provides: -- local path or `s3://` (MinIO), via `dissmodel.io.raster.save_geotiff` -- then `record.metrics.update(result["metrics"])`, output checksum (of the written file), status, final log. `crs`/`transform` come - straight off the `RasterBackend` (`load()` already sets both via - `vector_to_raster_backend`). Same mechanism as + straight off the `RasterBackend` (`load()` sets both, from the + GeoTIFF or via `vector_to_raster_backend`). Same mechanism as `brmangue-dissmodel`'s `RasterExecutor.save()`. `run(data, record)` remains abstract -- `data` arrives as the @@ -343,6 +439,29 @@ run the model on top of it and return `static`/`complementar_lu`/`allocation_data` for `transition_matrix` (`list[list[list[int]]]`, `[region][from][to]`). +### `LuccSaturationExecutor` + +`name = "lucc_continuous_saturation"`: `DemandPreComputedValues` + +`PotentialSpatialLagRegression` + `AllocationClueLikeSaturation`, the +continuous model of LuccME-BR. Parameters as the components', plus: + +- each `potential_data`/`allocation_data` entry names its class (`lu`) and, + optionally, its `region` (default 1) — one entry per class per region; +- `record.source.uri` is usually a **GeoTIFF with named bands** (read by the + base's `load()`): land uses, drivers and, optionally, `mask`, + `region`/`regionAloc` and the cell order (`order_attr`); +- `save_steps: list[int]` writes those steps too, as `_step.tif`. + +`examples/dissmodel-configs/lucc_saturation.toml` is LuccME's lab03; +`tests/test_executor_saturation.py` runs it through the local CLI and matches +the lab03 golden: + +```bash +python -m disslucc.executors.saturation run \ + --toml examples/dissmodel-configs/lucc_saturation.toml \ + --input cellspace.tif --param demand_csv=demand.csv --output out.tif +``` + ### Example ```python diff --git a/docs/architecture.md b/docs/architecture.md index efdefef..4e14c8a 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -123,6 +123,8 @@ scripts -- identical result (same checksum, same MAE, same F1). | `components/allocation/clue.py` | Yes, identical algorithm (main) | `disslucc-continuous/components/allocation/raster/clue.py` | | `components/potential/logistic.py` | Faithful algorithm, **raster is new** | ported from `disslucc-discrete/components/potential/vector/logistic_regression.py` (only vector existed) | | `components/allocation/clue_s.py` | Faithful algorithm, **raster is new** | ported from `disslucc-discrete/components/allocation/vector/clue_s.py` (only vector existed) | +| `components/potential/spatial_lag.py` | Yes, the Lua itself (LuccME `6244dd4`), **raster is new** | `luccme/lua/PotentialCSpatialLagRegression.lua` — no earlier Python port | +| `components/allocation/saturation.py` | Yes, the Lua itself, quirks included (`BACKP` across cells), **raster is new** | `luccme/lua/AllocationCClueLikeSaturation.lua` — no earlier Python port | | `validation/pontius.py` | New generalization | inspired by `disslucc-discrete/executors/lucc_validation_executor.py` (only had a binarized version) and matches the formula in `disslucc-continuous/executors/lucc_benchmark_executor.py::_metrics` (continuous, also not ported before) | **Detail that only shows up when comparing `disslucc-continuous` diff --git a/docs/decisions.md b/docs/decisions.md index d78a75a..6310122 100644 --- a/docs/decisions.md +++ b/docs/decisions.md @@ -1018,3 +1018,86 @@ reproduces TerraME year by year, iteration counts included (0, 0, 8, 26, 18, 17, 17; MAE < 1e-7). The tests against the continuous goldens use `False`; `test_lab1_default_cell_correction_deviates_from_terrame` pins the default's deviation. Both allocations also record `iterations_per_step`. + +## Spatial-lag potential and saturation allocation, for LuccME-BR (2026-09-24, #6) + +**What and why.** `PotentialSpatialLagRegression` and +`AllocationClueLikeSaturation` port LuccME's `PotentialCSpatialLagRegression` +and `AllocationCClueLikeSaturation` (LuccME `6244dd4`), the components LuccME-BR +(Bezerra et al. 2022, PLOS ONE e0256052) is built from. They were written and +validated first in [profsergiocosta/luccmebr-reconstruction](https://github.com/profsergiocosta/luccmebr-reconstruction), +with this package's interfaces, and moved here so that repository keeps only +the model. Nothing existing changed: Lab1/Lab15 numbers are the same. + +**Faithful to the Lua, including what looks accidental** — each point is pinned +by a test that fails on the other reading: +- the potential's constant adaptation accumulates year to year (LuccME writes + it back into `const`), unlike `PotentialLinearRegression`; +- `lab06`'s `updateYears` copies the new drivers into the cells by position + (`forEachCellPair`), and `csAC_2009` is not in `csAC`'s order; +- `correctCellChange` runs here (the Saturation variant spells `regionAloc` + right), as the Lua writes it — not `AllocationClueLike`'s own version — + and its `BACKP`, declared outside the loop over cells, carries over from cell + to cell: the cells' order matters, given by `order_attr`; +- the neighbourhood LuccME names "11x11" is `createNeighborhood{strategy="mxn"}` + with no `m`/`n`: TerraME's default, 3 × 3 with the cell + (`packages/base/lua/CellularSpace.lua`). + +**Validation, and the Lua in `tests/lua/`.** The `lab03`/`lab06` goldens match +cell for cell, year by year, iterations and maximum error included +(`docs/validation.md`). But those labs never reach `correctCellChange`, the +saturation branch or an isolated cell. For those, `tests/test_lua_differential.py` +runs the original Lua functions with `lupa` on synthetic cases that reach every +branch. That brings Lua source back into this repository, after the reference +scripts left it for terrame-docker (entry above) — deliberately, and +**provisionally**: these are component sources needed to test branches, not +reference scripts, and `lupa` + stubs is not TerraME. The intended replacement +is component-level goldens generated in terrame-docker (the same synthetic +cases run in the real TerraME, inputs and outputs as CSV); when they exist, +`tests/lua/` and the `lupa` dependency go, and `benchmark/goldens/` holds +results only again. + +**Left open.** `AllocationClueLike`'s `cell_correction=True` is an intended +version of the step LuccME skips there; the Saturation variant's +`correctCellChange` is LuccME's own. Whether the first should follow the second +is a separate question, not settled here. + +## `LuccSaturationExecutor`, with raster input (2026-09-24, #6) + +LuccME-BR is run from a model TOML, not a script: the same second entry point +the continuous and discrete executors are (entry "The Executor came back"). +`LuccSaturationExecutor` builds the three components from `record.parameters` +and differs from `LuccContinuousExecutor` where that model needs it: + +- **raster input.** A continental cellular space (≈ 260 k cells, 15+ drivers) + is built once, as a raster, by the data pipeline (DisSCube in + luccmebr-reconstruction); rasterizing a vector at every run, as the base's + `load()` does, is the wrong way round there. So `load()` reads a GeoTIFF + whose bands carry their names (`dissmodel.io.load_dataset(fmt="raster")`, + the format `save()` already writes), and falls back to the base for vectors. + `_read_geotiff` returns the georeference in `meta`, not in the backend, so + `load()` copies it over — `save()` needs it. +- **regions** are named per entry (`lu`, `region`) instead of implied by list + position, because LuccME-BR has three regions × six classes. +- **`save_steps`**: a model checked against maps of intermediate years + (LuccME-BR: IBGE 2010, 2012, 2014) needs those states, not only the last. + +Checked end to end through the CLI: `examples/dissmodel-configs/lucc_saturation.toml` +is lab03, and its run matches the lab03 golden in 2011 and 2014 +(`tests/test_executor_saturation.py`). + +## GeoTIFF input moved to the base executor (2026-09-24, #6) + +The named-band GeoTIFF reading `LuccSaturationExecutor.load()` introduced +(entry above) is not specific to that model: any LUCC executor whose +cellular space is built once by a data pipeline wants it. It now lives in +`LuccExecutorBase.load()` -- a `.tif`/`.tiff` source is read as it is, +anything else is rasterized as before -- and `LuccSaturationExecutor` no +longer overrides `load()`. The georeference is still copied from `meta`: +dissmodel releases up to 0.6.5 keep it only there. + +In the same pass, the spatial lag's private shift helper gave way to +`dissmodel`'s `RasterBackend.shift2d`: same slices, opposite sign +convention, and the Moore neighbourhood is symmetric, so the sums are the +same (the golden and Lua-differential tests pass unchanged). + diff --git a/docs/validation.md b/docs/validation.md index 132da50..152c497 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -121,9 +121,9 @@ against the TerraME log -- done now, year by year, in the next section. for every cell and every simulated year, `_out` and `_pot`, plus TerraME's convergence-loop iteration count per year. terrame-docker has 23 (the 21 functional labs of the LuccME package and the two scenarios above); -this repository keeps the four its tests use (`lab01`, `lab01_md1643`, -`lab15`, `lab15_md10`), and adds the others as their components are -implemented. +this repository keeps the ones its tests use (`lab01`, `lab01_md1643`, +`lab15`, `lab15_md10`, and `lab03`/`lab06` for the components of the next +section), and adds the others as their components are implemented. `tests/test_goldens_per_year.py` checks the iteration count per year (exact) and every class per year (MAE < 1e-6): @@ -142,6 +142,35 @@ This closes the gap left by the Lab15 discriminance warning: the final map alone could not tell a correct CLUE-S from a static ranking, but the iteration counts can, and they match in every year. +## Spatial-lag potential and saturation allocation (lab03, lab06) + +`PotentialSpatialLagRegression` and `AllocationClueLikeSaturation` port the two +LuccME components LuccME-BR (Bezerra et al. 2022) is built from. `lab03` and +`lab06` are the two LuccME labs that combine exactly +`DemandPreComputedValues` + `PotentialCSpatialLagRegression` + +`AllocationCClueLikeSaturation` (csAC, 2008–2014); `lab06` also replaces the +driver `ti` in 2009 (`updateYears`). + +| Check | Result | Test | +|---|---|---| +| potential alone, each year from the golden's previous land use | every cell, year and class, \|Δ\| < 1e-10 (the goldens keep 12 decimals) | `test_spatial_lag_golden.py` | +| whole model (potential + allocation) run from 2008 on its own | every cell's `_out` and `_pot`, every year to 2014, \|Δ\| < 1e-9 | `test_saturation_golden.py` | +| TerraME log, per year | iterations (0 every year) and "Maximum error" identical (rel. 1e-9) | `test_saturation_golden.py` | + +What the goldens also tell apart: a non-cumulative constant adaptation (the +reading `PotentialLinearRegression` uses) misses from 2010 on; `lab06`'s +dynamic variables paired by `object_id0` instead of by position (LuccME's +`forEachCellPair`; `csAC_2009` is not in `csAC`'s order) misses by 0.045. + +What they do **not** reach: no cell ever needs `correctCellChange`, no cell is +saturated (`changeLimiarValue = 1`), and csAC has no isolated cell — the cells' +visiting order has no effect either. Those parts are checked against the +original Lua functions (`tests/lua/`, run with `lupa`) on synthetic cases that +reach every branch, within 1e-12 (`test_lua_differential.py`); ports broken on +purpose (`BACKP` reset per cell, the 3 × 3 window without the cell, an isolated +cell using its own share) fail them. Not tested at all: `modify_driver`, called +only after 500 iterations. + ## Summary | Scenario | Type | MAE | Match | What it proves | @@ -149,6 +178,7 @@ iteration counts can, and they match in every year. | Lab1 | continuous | 0.0036 | -- | full model (Demand+Potential+Allocation), within the official tolerance; the whole MAE is the cell correction that TerraME skips | | Lab15 | discrete | 0.0 | 100% | correct regression coefficients | | Lab1, Lab15 year by year | both | < 1e-7 | iterations exact | CLUE-S convergence confirmed; CLUE identical to TerraME with `cell_correction=False` | +| lab03, lab06 year by year | continuous | < 1e-9 | iterations and max error exact | spatial-lag potential + saturation allocation identical to TerraME; the branches the labs don't reach, identical to the Lua | ## What these numbers DON'T prove: engineering validation ≠ scientific validation diff --git a/examples/dissmodel-configs/lucc_saturation.toml b/examples/dissmodel-configs/lucc_saturation.toml new file mode 100644 index 0000000..e13e540 --- /dev/null +++ b/examples/dissmodel-configs/lucc_saturation.toml @@ -0,0 +1,97 @@ +# configs/models/lucc_saturation.toml +# +# Registers disslucc's LuccSaturationExecutor (name="lucc_continuous_saturation"): +# DemandPreComputedValues + PotentialSpatialLagRegression + +# AllocationClueLikeSaturation, the components of LuccME-BR. +# +# This file is LuccME's lab03 (luccme/tests/functional/lab03.lua), the same +# coefficients, demand and allocation data; tests/test_executor_saturation.py +# runs it through the local CLI and checks the result against the lab03 golden: +# +# python -m disslucc.executors.saturation run \ +# --toml examples/dissmodel-configs/lucc_saturation.toml \ +# --input \ +# --param demand_csv=data/input/examples_demand_lab1.csv --output out.tif +# +# Every [[model.potential_data]] / [[model.allocation_data]] entry names its +# class (lu) and may name its region (default 1): one entry per class per region. + +[model] +executor_module = "disslucc.executors" +name = "lucc_continuous_saturation" +class = "lucc_continuous_saturation" +description = "LUCC continuous simulation (CLUE-like with saturation, spatial-lag regression potential)" +package = "git+https://github.com/DisSModel/disslucc@main" +dissmodel = ">=0.6.4,<0.7.0" + +land_use_types = ["f", "d", "outros"] +land_use_no_data = "outros" +complementar_lu = "f" +cell_area = 25.0 +attr_protection = "uc_pi" +saturation_indicator = "saturationLimiar" +order_attr = "order" # LuccME visits cells in the layer's order (BACKP in correctCellChange) +max_difference = 1643.0 +max_iteration = 1000 +initial_elasticity = 0.1 +min_elasticity = 0.001 +max_elasticity = 1.5 + +# Runtime parameters. demand_csv comes with the experiment. +[model.parameters] +n_steps = 7 +save_steps = [3] # also write step 3 (2011) + +[[model.potential_data]] +lu = "f" +const = 0.05266679 +ro = 0.9124615 + [model.potential_data.betas] + uc_us = 0.03789872 + uc_pi = 0.04141921 + ti = 0.04455667 + +[[model.potential_data]] +lu = "d" +const = 0.01431553 +ro = 0.9019253 + [model.potential_data.betas] + assentamen = 0.0443537 + uc_us = -0.01454847 + fertilidad = 0.01701601 + dist_riobr = -0.00000002262071 + +[[model.potential_data]] +lu = "outros" +const = 0.0 +ro = 0.0 + +[[model.allocation_data]] +lu = "f" +static = -1 +min_value = 0.0 +max_value = 1.0 +min_change = 0.0 +max_change = 1.0 +change_limiar_value = 1.0 +max_change_above_limiar = 0.0 + +[[model.allocation_data]] +lu = "d" +static = -1 +min_value = 0.0 +max_value = 1.0 +min_change = 0.0 +max_change = 1.0 +change_limiar_value = 1.0 +max_change_above_limiar = 0.0 + +[[model.allocation_data]] +lu = "outros" +static = 1 +min_value = 0.0 +max_value = 1.0 +min_change = 0.0 +max_change = 1.0 +change_limiar_value = 1.0 +max_change_above_limiar = 0.0 diff --git a/pyproject.toml b/pyproject.toml index 12a8b95..f526f8a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,7 +17,9 @@ dependencies = [ [project.optional-dependencies] examples = ["matplotlib"] -dev = ["mypy>=1.10", "pytest>=7", "ruff>=0.16,<0.17"] +# lupa runs the original LuccME Lua in tests/test_lua_differential.py +# (skipped without it); see docs/decisions.md +dev = ["mypy>=1.10", "pytest>=7", "ruff>=0.16,<0.17", "lupa"] [tool.mypy] # 3.12, not 3.11 (the floor in requires-python): numpy 2.5's stub uses diff --git a/src/disslucc/__init__.py b/src/disslucc/__init__.py index d66b320..07d30ee 100644 --- a/src/disslucc/__init__.py +++ b/src/disslucc/__init__.py @@ -3,14 +3,21 @@ (main branch). See README for what was simplified and what wasn't (the CLUE algorithm itself is faithful to the original). """ -from .components.allocation import AllocationClueLike, AllocationDClueSLike +from .components.allocation import AllocationClueLike, AllocationClueLikeSaturation, AllocationDClueSLike from .components.demand import DemandInline, DemandPreComputedValues, load_demand_csv from .components.potential import ( PotentialDLogisticRegression, PotentialLinearRegression, + PotentialSpatialLagRegression, +) +from .protocols import DemandProtocol, PotentialProtocol, RegionalPotentialProtocol +from .schemas import ( + AllocationSpec, + LogisticRegressionSpec, + RegressionSpec, + SaturationAllocationSpec, + SpatialLagRegressionSpec, ) -from .protocols import DemandProtocol, PotentialProtocol -from .schemas import AllocationSpec, LogisticRegressionSpec, RegressionSpec from .validation import confusion_metrics, pontius_millones # Explicit public API: these are re-exports (nothing here is used within @@ -18,6 +25,7 @@ # every one of them as an unused import. __all__ = [ "AllocationClueLike", + "AllocationClueLikeSaturation", "AllocationDClueSLike", "AllocationSpec", "DemandInline", @@ -27,7 +35,11 @@ "PotentialDLogisticRegression", "PotentialLinearRegression", "PotentialProtocol", + "PotentialSpatialLagRegression", + "RegionalPotentialProtocol", "RegressionSpec", + "SaturationAllocationSpec", + "SpatialLagRegressionSpec", "confusion_metrics", "load_demand_csv", "pontius_millones", diff --git a/src/disslucc/components/allocation/__init__.py b/src/disslucc/components/allocation/__init__.py index f4f8ada..6e7c983 100644 --- a/src/disslucc/components/allocation/__init__.py +++ b/src/disslucc/components/allocation/__init__.py @@ -1,4 +1,5 @@ from .clue import AllocationClueLike from .clue_s import AllocationDClueSLike +from .saturation import AllocationClueLikeSaturation -__all__ = ["AllocationClueLike", "AllocationDClueSLike"] +__all__ = ["AllocationClueLike", "AllocationClueLikeSaturation", "AllocationDClueSLike"] diff --git a/src/disslucc/components/allocation/saturation.py b/src/disslucc/components/allocation/saturation.py new file mode 100644 index 0000000..f45095d --- /dev/null +++ b/src/disslucc/components/allocation/saturation.py @@ -0,0 +1,417 @@ +""" +disslucc.components.allocation.saturation +----------------------------------------- +Port of LuccME's ``AllocationCClueLikeSaturation`` (TerraME, LuccME 3.1, commit +6244dd4 — ``luccme/lua/AllocationCClueLikeSaturation.lua``), raster only. Used +by LuccME-BR (Bezerra et al. 2022). + +CLUE-like continuous allocation (Verburg et al. 1999): each year, the change of a +class in a cell is ``potential × elasticity``, bounded per class and region; the +elasticity of each class is adapted until the allocated areas meet the demand. +What the Saturation variant adds to ``AllocationCClueLike``: + +- a **saturation indicator** per cell, recomputed at the start of every year: the + share of the "available" area (not no-data, not protected) that is no longer in + the complementary class, averaged over the cell's 3 × 3 window. Where it + exceeds ``change_limiar_value``, the change in the demand's direction is halved + or capped at ``max_change_above_limiar``; +- ``correctCellChange`` **runs** (in ``AllocationCClueLike`` a typo in its region + test — ``cell.regionregionAloc`` — skips it, which is why + ``AllocationClueLike`` implements its own version, behind ``cell_correction``; + this one is the Lua's own ``correctCellChange``); +- small differences in the bounds (no cap at 1, ``< min_value`` instead of + ``<=``), a guard for a zero potential, areas counted over positive values only, + and ``potential.modify_driver`` when half the iterations have gone by. + +Faithful to the Lua, including what looks accidental: + +- ``BACKP`` in ``correctCellChange`` is declared outside the loop over cells: a + cell that lowers it lowers it for the cells visited after it, within the same + call. The result depends on the order of the cells — ``order_attr`` gives it + (TerraME visits the cells in the order of the layer); +- the neighbourhood LuccME names "11x11" is created with ``strategy = "mxn"`` and + no ``m``/``n``: TerraME's default, 3 × 3 including the cell itself; +- ``compareAllocationToDemand`` adapts the elasticities once per potential region, + and reads ``static`` from region 1. + +Validation: +- tests/test_saturation_golden.py — the whole model (this allocation + the + spatial-lag potential) run from 2008 against the TerraME goldens of lab03 and + lab06: every cell's land use and potential, every year to 2014, within 1e-9, + plus TerraME's logged iterations and maximum error per year. Those labs never + need ``correctCellChange`` and never reach the saturation branch; +- tests/test_lua_differential.py — ``correct_cell_change``, ``compute_change`` + (saturation branch included) and ``saturation_indicator`` against the + original Lua functions, run with lupa on synthetic cases that reach every + branch, the running ``BACKP`` included (within 1e-12). +Not tested: ``potential.modify_driver`` (called only after 500 iterations). +""" + +from __future__ import annotations + +from typing import cast + +import numpy as np +from dissmodel.geo import SyncRasterModel + +from ...protocols import DemandProtocol, RegionalPotentialProtocol +from ...schemas import SaturationAllocationSpec + +TOL_COVER = 0.005 # LuccME: a cell whose classes sum to 1 ± 0.005 is left alone +MAX_CORRECTIONS = 25 # LuccME: iterations of the per-cell correction +BACKP_START = 0.5 + + +def compute_change( + past: np.ndarray, + pot: np.ndarray, + elasticity: float, + spec: SaturationAllocationSpec, + direction: int, + saturation: np.ndarray, +) -> np.ndarray: + """New share of one class (LuccME ``computeChange``), elementwise.""" + change = pot * elasticity + small = np.abs(change) < spec.min_change + pot = np.where(small, 0.0, pot) + change = np.where(small, 0.0, change) + big = (np.abs(change) >= spec.max_change) & (pot != 0) + change = np.where(big, spec.max_change * np.sign(pot), change) + + saturated = saturation > spec.change_limiar_value + if spec.static < 1: + cap = spec.max_change_above_limiar + up = saturated & (pot >= 0) & (direction == 1) & (change >= cap) + change = np.where(up, np.where(change / 2 < cap, change / 2, cap), change) + down = saturated & (pot <= 0) & (direction == -1) & (np.abs(change) >= cap) + change = np.where(down, np.where(np.abs(change / 2) < cap, change / 2, -cap), change) + + if spec.static == 1: + new = past.copy() + elif spec.static == 0: + new = past + change + else: + follows = ((pot >= 0) & (direction == 1)) | ((pot <= 0) & (direction == -1)) + new = np.where(follows, past + change, past) + + new = np.where(new < 0, 0.0, new) + new = np.where(new < spec.min_value, np.where(past >= spec.min_value, spec.min_value, past), new) + new = np.where(new > spec.max_value, np.where(past <= spec.max_value, spec.max_value, past), new) + return new + + +def correct_cell_change( + values: np.ndarray, + past: np.ndarray, + specs: list[SaturationAllocationSpec], + directions: list[int], + stats: dict | None = None, +) -> np.ndarray: + """LuccME ``correctCellChange``: bring each cell's classes back to a total of 1. + + ``values``, ``past``: (cells, classes), the cells **in visiting order** — the + running ``BACKP`` depends on it. Returns the corrected values. ``stats``, if + given, receives how many cells took each branch (for tests). + """ + v = values.astype(np.float64).copy() + p = past.astype(np.float64) + static = np.array([s.static for s in specs]) + vmin = np.array([s.min_value for s in specs]) + vmax = np.array([s.max_value for s in specs]) + direction = np.array(directions) + + need = np.abs(v.sum(axis=1) - 1) > TOL_COVER + dif = v - p + totchange = np.abs(dif).sum(axis=1) + biggest = np.abs(dif).max(axis=1) + amin, amax = dif.min(axis=1), dif.max(axis=1) + + # BACKP carries over from cell to cell: a running minimum in visiting order + with np.errstate(divide="ignore", invalid="ignore"): + candidate = np.where(need & (totchange > 0), biggest / (2 * totchange), np.inf) + backp = np.minimum(BACKP_START, np.minimum.accumulate(candidate)) + + def free(x: np.ndarray) -> np.ndarray: + return (static < 1) & ~((x <= vmin) | (x >= vmax)) + + # all free classes moving the same way: shift them back + f = free(v) + step = (backp * totchange)[:, None] + incr = (f & (dif > -step)).sum(axis=1) + decr = (f & (dif < step)).sum(axis=1) + n_free = f.sum(axis=1) + all_incr = need & (incr == n_free) + all_decr = need & (decr == n_free) + shift = all_incr | all_decr + if shift.any(): + s = v.copy() + s = np.where(f & all_incr[:, None], s - (amin + backp * totchange)[:, None], s) + s = np.where(f & all_decr[:, None], s + (backp * totchange - amax)[:, None], s) + s = np.where(f & (s < 0), 0.0, s) + follows = static == -1 + s = np.where(follows & (direction == 1) & (s < p), p, s) + s = np.where(follows & (direction == -1) & (s > p), p, s) + v = np.where(shift[:, None], s, v) + + # proportional correction, up to 25 times per cell + active = need.copy() + rounds = np.zeros(len(v), dtype=np.int64) + for _ in range(MAX_CORRECTIONS): + if not active.any(): + break + rounds += active + f = free(v) + totcov = (v * f).sum(axis=1) + totstatic = (v * ~f).sum(axis=1) + totch = np.abs(v - p).sum(axis=1) + off = active & (np.abs(totcov - (1 - totstatic)) > TOL_COVER) + with np.errstate(divide="ignore", invalid="ignore"): + rescale = np.where(f, v * ((1 - totstatic) / totcov)[:, None], v) + reduce = v - np.abs(v - p) * ((totcov - (1 - totstatic)) / totch)[:, None] + reduce = np.where(reduce < 0, 0.0, reduce) + v = np.where((off & (totch == 0))[:, None], rescale, v) + v = np.where((off & (totch != 0))[:, None], reduce, v) + active = off # the Lua stops a cell once its check passes at the start of a round + + last = need & (rounds == MAX_CORRECTIONS) + if stats is not None: + stats["need"] = stats.get("need", 0) + int(need.sum()) + stats["shift"] = stats.get("shift", 0) + int(shift.sum()) + stats["backp_lowered"] = stats.get("backp_lowered", 0) + int((need & (backp < BACKP_START)).sum()) + stats["several_rounds"] = stats.get("several_rounds", 0) + int((rounds > 1).sum()) + stats["last"] = stats.get("last", 0) + int(last.sum()) + if last.any(): + f = free(v) + totcov = (v * f).sum(axis=1) + totstatic = (v * ~f).sum(axis=1) + with np.errstate(divide="ignore", invalid="ignore"): + v = np.where(last[:, None] & f, v * ((1 - totstatic) / totcov)[:, None], v) + return v + + +def saturation_indicator( + complementar: np.ndarray, + valid: np.ndarray, + no_data: np.ndarray | None = None, + protection: np.ndarray | None = None, +) -> np.ndarray: + """LuccME ``updateAllocationParameters``: share of the available area no longer + in the complementary class, averaged over the 3 × 3 window (cell included).""" + original = np.ones_like(complementar) if no_data is None else 1 - no_data + prot = np.zeros_like(complementar) if protection is None else protection + available = original - prot + with np.errstate(divide="ignore", invalid="ignore"): + perc = np.minimum((1 - complementar) / available, 1.0) + ok = valid & (available > 0) + + total = np.zeros_like(complementar, dtype=np.float64) + count = np.zeros(complementar.shape, dtype=np.int64) + rows, cols = complementar.shape + for dr in (-1, 0, 1): + for dc in (-1, 0, 1): + r0, r1 = max(0, -dr), min(rows, rows - dr) + c0, c1 = max(0, -dc), min(cols, cols - dc) + src_ok = ok[r0 + dr : r1 + dr, c0 + dc : c1 + dc] + total[r0:r1, c0:c1] += np.where(src_ok, perc[r0 + dr : r1 + dr, c0 + dc : c1 + dc], 0.0) + count[r0:r1, c0:c1] += src_ok + with np.errstate(divide="ignore", invalid="ignore"): + averaged = np.where(count > 0, total / count, perc) + return np.where(available > 0, averaged, 1.0) + + +class AllocationClueLikeSaturation(SyncRasterModel): + """Continuous CLUE allocation with saturation — interface of ``AllocationClueLike``. + + Arrays in the backend: the land uses (their ``_past`` kept by + ``SyncRasterModel``), ``_pot`` written by the potential, ``mask`` (cells + that exist), ``region_attr`` (default: all cells in region 1), and optionally + ``order_attr`` (the order in which LuccME visits the cells; default: row-major) + and ``attr_protection``. Writes the land uses, ``_out`` and the saturation + indicator. + """ + + def setup( + self, + backend, + demand: DemandProtocol, + potential: RegionalPotentialProtocol, + land_use_types: list[str], + allocation_data: list[list[SaturationAllocationSpec]], + complementar_lu: str, + cell_area: float, + land_use_no_data: str | None = None, + attr_protection: str | None = None, + saturation_indicator: str = "saturationLimiar", + max_difference: float = 1643, + max_iteration: int = 1000, + initial_elasticity: float = 0.1, + min_elasticity: float = 0.001, + max_elasticity: float = 1.5, + region_attr: str = "regionAloc", + mask_attr: str = "mask", + order_attr: str | None = None, + ) -> None: + super().setup(backend) + self.demand = demand + self.potential = potential + self.land_use_types = land_use_types + self.allocation_data = allocation_data + self.complementar_lu = complementar_lu + self.cell_area = cell_area + self.land_use_no_data = land_use_no_data + self.attr_protection = attr_protection + self.saturation_attr = saturation_indicator + self.max_difference = max_difference + self.max_iteration = max_iteration + self.initial_elasticity = initial_elasticity + self.min_elasticity = min_elasticity + self.max_elasticity = max_elasticity + self.region_attr = region_attr + self.mask_attr = mask_attr + self.order_attr = order_attr + # per step, as LuccME logs them: iterations of the convergence loop + # (0 = first allocation accepted) and the maximum |area - demand| accepted + self.iterations_per_step: list[int] = [] + self.max_error_per_step: list[float] = [] + + if region_attr not in self.backend.arrays: + self.backend.set(region_attr, np.ones(self.shape, dtype=np.int32)) + if mask_attr not in self.backend.arrays: + self.backend.set(mask_attr, np.ones(self.shape, dtype=np.float32)) + + # ── helpers ───────────────────────────────────────────────────────────── + + def _valid(self) -> np.ndarray: + return cast(np.ndarray, self.backend.get(self.mask_attr) > 0) + + def _cells(self, r_number: int) -> tuple[np.ndarray, np.ndarray]: + """(rows, cols) of the region's cells, in LuccME's visiting order.""" + sel = self._valid() & (self.backend.get(self.region_attr) == r_number) + rows, cols = np.nonzero(sel) + if self.order_attr is not None: + key = self.backend.get(self.order_attr)[rows, cols] + idx = np.argsort(key, kind="stable") + rows, cols = rows[idx], cols[idx] + return rows, cols + + def _areas(self) -> list[float]: + valid = self._valid() + return [ + float(np.where(self.backend.get(lu) > 0, self.backend.get(lu), 0.0)[valid].sum()) * self.cell_area + for lu in self.land_use_types + ] + + # ── the year ──────────────────────────────────────────────────────────── + + def execute(self) -> None: + step = int(self.env.now()) + self.elasticity = [self.initial_elasticity] * len(self.land_use_types) + self.update_saturation() + + n_iter, max_adjust, flex = 0, self.max_difference, False + while True: + if step != 0: + for r_number in range(1, len(self.allocation_data) + 1): + self.compute_change(r_number) + self.correct_cell_change(r_number) + max_diff = self.compare_to_demand(step) + if max_diff <= max_adjust: + break + n_iter += 1 + if n_iter > self.max_iteration * 0.5 and not flex: + max_adjust *= 2 + flex = True + if self.attr_protection is not None: + self.potential.modify_driver(self.attr_protection, 0.5) + if n_iter >= self.max_iteration: + raise RuntimeError(f"allocation did not converge at step {step} (error {max_diff:.1f})") + self.iterations_per_step.append(n_iter) + self.max_error_per_step.append(max_diff) + + self.apply_complementar() + for lu in self.land_use_types: + self.backend.set(lu + "_out", self.backend.get(lu).copy()) + + def update_saturation(self) -> None: + get = self.backend.get + self.backend.set( + self.saturation_attr, + saturation_indicator( + get(self.complementar_lu).astype(np.float64), + self._valid(), + get(self.land_use_no_data).astype(np.float64) if self.land_use_no_data else None, + get(self.attr_protection).astype(np.float64) if self.attr_protection else None, + ), + ) + + def compute_change(self, r_number: int) -> None: + sel = self._valid() & (self.backend.get(self.region_attr) == r_number) + saturation = self.backend.get(self.saturation_attr) + for lu_idx, lu in enumerate(self.land_use_types): + spec = self.allocation_data[r_number - 1][lu_idx] + new = compute_change( + self.backend.get(lu + "_past").astype(np.float64), + self.backend.get(lu + "_pot").astype(np.float64), + self.elasticity[lu_idx], + spec, + self.demand.get_current_lu_direction(lu_idx), + saturation, + ) + self.backend.set(lu, np.where(sel, new, self.backend.get(lu))) + + def correct_cell_change(self, r_number: int) -> None: + rows, cols = self._cells(r_number) + lus = self.land_use_types + values = np.stack([self.backend.get(lu)[rows, cols] for lu in lus], axis=1) + past = np.stack([self.backend.get(lu + "_past")[rows, cols] for lu in lus], axis=1) + directions = [self.demand.get_current_lu_direction(i) for i in range(len(lus))] + corrected = correct_cell_change(values, past, self.allocation_data[r_number - 1], directions) + for i, lu in enumerate(lus): + arr = self.backend.get(lu).astype(np.float64).copy() + arr[rows, cols] = corrected[:, i] + self.backend.set(lu, arr) + + def compare_to_demand(self, step: int) -> float: + areas = self._areas() + max_diff = 0.0 + for j in range(1, len(self.potential.potential_data) + 1): + for i in range(len(self.land_use_types)): + direction = self.demand.get_current_lu_direction(i) + demand = self.demand.get_current_lu_demand(i) + if direction == 0 and step == 0: + direction = 1 if demand >= areas[i] else -1 + if direction == 1: + self.elasticity[i] *= demand / areas[i] + else: + self.elasticity[i] *= areas[i] / demand + if self.elasticity[i] > self.max_elasticity: + self.elasticity[i] = self.max_elasticity + self.potential.modify(j, i, direction) + if self.elasticity[i] < self.min_elasticity: + if self.allocation_data[0][i].static < 0: + self.elasticity[i] = self.min_elasticity + self.potential.modify(j, i, -direction) + else: + self.demand.change_lu_direction(i) + max_diff = max(max_diff, abs(areas[i] - demand)) + return max_diff + + def apply_complementar(self) -> None: + """complementar = 1 − the others; a deficit comes off the largest class + (neither the complementary nor the no-data one).""" + valid = self._valid() + others = [lu for lu in self.land_use_types if lu != self.complementar_lu] + total = sum(self.backend.get(lu).astype(np.float64) for lu in others) + comp = 1 - total + + candidates = [lu for lu in others if lu != self.land_use_no_data] + if candidates: + stack = np.stack([self.backend.get(lu).astype(np.float64) for lu in candidates]) + biggest = np.argmax(stack, axis=0) # first maximum, as LuccME's strict ">" + has_bigger = stack.max(axis=0) > 0 + deficit = valid & (comp < 0) & has_bigger + for k, lu in enumerate(candidates): + hit = deficit & (biggest == k) + self.backend.set(lu, np.where(hit, self.backend.get(lu) + comp, self.backend.get(lu))) + comp = np.where(comp < 0, 0.0, comp) + self.backend.set(self.complementar_lu, np.where(valid, comp, self.backend.get(self.complementar_lu))) diff --git a/src/disslucc/components/potential/__init__.py b/src/disslucc/components/potential/__init__.py index 930e7c4..b4c7bf3 100644 --- a/src/disslucc/components/potential/__init__.py +++ b/src/disslucc/components/potential/__init__.py @@ -1,4 +1,5 @@ from .linear import PotentialLinearRegression from .logistic import PotentialDLogisticRegression +from .spatial_lag import PotentialSpatialLagRegression -__all__ = ["PotentialDLogisticRegression", "PotentialLinearRegression"] +__all__ = ["PotentialDLogisticRegression", "PotentialLinearRegression", "PotentialSpatialLagRegression"] diff --git a/src/disslucc/components/potential/spatial_lag.py b/src/disslucc/components/potential/spatial_lag.py new file mode 100644 index 0000000..2eed9d9 --- /dev/null +++ b/src/disslucc/components/potential/spatial_lag.py @@ -0,0 +1,203 @@ +""" +disslucc.components.potential.spatial_lag +----------------------------------------- +Port of LuccME's ``PotentialCSpatialLagRegression`` (TerraME, LuccME 3.1, +commit 6244dd4 — ``luccme/lua/PotentialCSpatialLagRegression.lua``), raster +only. Used by LuccME-BR (Bezerra et al. 2022) for every class. + +For each class, in every cell of a region:: + + Y = mean of the class share over the cell and its Moore neighbours, + each share divided by (1 - no-data share); neighbours that are + all no-data left out + reg = newconst + Σ beta·x + ro·Y (log10 space if is_log) + limit = const + Σ beta·x + ro·Y + reg = 1 if limit > max_reg, 0 if limit < min_reg + reg = reg · (1 - no-data share) + pot = reg - past share + +Quirks of the original, kept on purpose: +- the clip tests ``limit`` (built with ``const``) but writes 1 or 0, not + ``max_reg`` / ``min_reg`` — exercised by the goldens; +- ``adapt_constants`` writes the adapted constant back into ``const``, so the + adaptation accumulates year after year — the goldens reject the other reading; +- a cell with no neighbour at all gets Y = 0, not its own share — no cell of + csAC is isolated; checked against the Lua itself instead. + +Validated against the LuccME goldens of lab03 and lab06 (2008–2014, every +cell, |Δ| < 1e-10: tests/test_spatial_lag_golden.py) and against the original +Lua function run with lupa, log-transformed classes, isolated and all-no-data +cells included (tests/test_lua_differential.py). +""" + +from __future__ import annotations + +import math + +import numpy as np +from dissmodel.geo import RasterBackend, SyncRasterModel + +from ...protocols import DemandProtocol +from ...schemas import SpatialLagRegressionSpec + +MOORE = [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)] +LOG_OFFSET = 0.0001 # LuccME's "ANAP" offset for log-transformed classes +CONST_CHANGE = 0.1 # LuccME's constChange: the allocation's step on newconst + + +def spatial_lag_regression( + spec: SpatialLagRegressionSpec, + past: np.ndarray, + drivers: dict[str, np.ndarray], + valid: np.ndarray, + no_data: np.ndarray | None = None, + newconst: float | None = None, +) -> tuple[np.ndarray, np.ndarray]: + """``(reg, pot)`` of one class, on every cell (callers apply the region mask). + + ``past``: share of the class at the start of the step; ``valid``: cells + that exist (LuccME cells); ``no_data``: share of the no-data class, or + None; ``newconst``: defaults to ``spec.const``. + """ + past = np.asarray(past, dtype=np.float64) + valid = np.asarray(valid, dtype=bool) + nd = np.zeros_like(past) if no_data is None else np.asarray(no_data, dtype=np.float64) + newconst = spec.const if newconst is None else newconst + + with np.errstate(divide="ignore", invalid="ignore"): + scaled = np.where(nd != 1, past / (1 - nd), past) + eligible = valid & (nd < 1) + + neigh_sum = np.zeros_like(past) + count = np.zeros(past.shape, dtype=np.int64) + exists = np.zeros(past.shape, dtype=np.int64) + for dr, dc in MOORE: + e = RasterBackend.shift2d(eligible, dr, dc) + neigh_sum += np.where(e, RasterBackend.shift2d(scaled, dr, dc), 0.0) + count += e + exists += RasterBackend.shift2d(valid, dr, dc) + + y = np.where(count > 0, (scaled + neigh_sum) / (count + 1), np.where(exists > 0, scaled, 0.0)) + if spec.is_log: + y = np.log10(y + LOG_OFFSET) + y = y * spec.ro + + x = np.zeros_like(past) + for name, beta in spec.betas.items(): + x = x + beta * np.asarray(drivers[name], dtype=np.float64) + + reg = newconst + x + y + limit = spec.const + x + y + if spec.is_log: + reg = 10.0**reg - LOG_OFFSET + limit = 10.0**limit - LOG_OFFSET + reg = np.where(limit > spec.max_reg, 1.0, reg) + reg = np.where(limit < spec.min_reg, 0.0, reg) + reg = reg * (1 - nd) + return reg, reg - past + + +class PotentialSpatialLagRegression(SyncRasterModel): + """Potential by spatial-lag regression — interface of ``PotentialLinearRegression``. + + Arrays in the backend: the land uses, their ``_past`` (kept by + ``SyncRasterModel``), the drivers named in the betas, ``region_attr`` + (default: all cells in region 1) and ``mask`` (cells that exist). + Writes ``_pot`` and ``_reg``. + """ + + def setup( + self, + backend, + potential_data: list[list[SpatialLagRegressionSpec]], + demand: DemandProtocol, + land_use_types: list[str], + land_use_no_data: str | None = None, + region_attr: str = "region", + mask_attr: str = "mask", + ) -> None: + super().setup(backend) + self.potential_data = potential_data + self.demand = demand + self.land_use_types = land_use_types + self.land_use_no_data = land_use_no_data + self.region_attr = region_attr + self.mask_attr = mask_attr + + if region_attr not in self.backend.arrays: + self.backend.set(region_attr, np.ones(self.shape, dtype=np.int32)) + if mask_attr not in self.backend.arrays: + self.backend.set(mask_attr, np.ones(self.shape, dtype=np.float32)) + for region in self.potential_data: + for spec in region: + spec.newconst = spec.const + for lu in self.land_use_types: + self.backend.set(lu + "_pot", np.zeros(self.shape, dtype=np.float64)) + self.backend.set(lu + "_reg", np.zeros(self.shape, dtype=np.float64)) + + def execute(self) -> None: + step = int(self.env.now()) + for r_idx, region in enumerate(self.potential_data): + for spec in region: + spec.newconst = spec.const + if step > 0: + self.adapt_constants(r_idx + 1) + for lu_idx in range(len(self.land_use_types)): + self.compute_potential(r_idx + 1, lu_idx) + + def adapt_constants(self, r_number: int) -> None: + """LuccME ``adaptRegressionConstants``: const += 0.01 · relative change of demand.""" + for lu_idx, spec in enumerate(self.potential_data[r_number - 1]): + curr = self.demand.get_current_lu_demand(lu_idx) + prev = self.demand.get_previous_lu_demand(lu_idx) + plus = 0.01 * ((curr - prev) / prev) + spec.newconst = spec.const + if spec.is_log: + unlog = 10**spec.newconst + plus + if unlog != 0: + spec.newconst = math.log10(unlog) + else: + spec.newconst = spec.newconst + plus + spec.const = spec.newconst + + def compute_potential(self, r_number: int, lu_idx: int) -> None: + lu = self.land_use_types[lu_idx] + spec = self.potential_data[r_number - 1][lu_idx] + in_region = self.backend.get(self.region_attr) == r_number + drivers = {name: self.backend.get(name) for name in spec.betas} + no_data = self.backend.get(self.land_use_no_data) if self.land_use_no_data else None + reg, pot = spatial_lag_regression( + spec, + self.backend.get(lu + "_past"), + drivers, + self.backend.get(self.mask_attr) > 0, + no_data, + spec.newconst, + ) + self.backend.arrays[lu + "_reg"] = np.where(in_region, reg, self.backend.get(lu + "_reg")) + self.backend.arrays[lu + "_pot"] = np.where(in_region, pot, self.backend.get(lu + "_pot")) + + def modify_driver(self, attr_protection: str, rate: float) -> None: + """LuccME ``modifyDriver``: scale the protection driver's beta in every regression. + + Called by the saturation allocation when it has not converged after half + its iterations. As in the Lua, only the first class is recomputed: the loop + meant to find the complementary class shadows its own variable and always + stops at index 1. Not exercised by the lab03/lab06 goldens. + """ + for r_idx, region in enumerate(self.potential_data): + for spec in region: + if attr_protection in spec.betas: + spec.betas[attr_protection] *= rate + self.compute_potential(r_idx + 1, 0) + + def modify(self, r_number: int, lu_idx: int, direction: int) -> None: + """LuccME ``modify``: the allocation moves newconst by ±constChange and recomputes.""" + spec = self.potential_data[r_number - 1][lu_idx] + if spec.is_log: + unlog = 10**spec.newconst + CONST_CHANGE * direction + if unlog != 0: + spec.newconst = math.log10(unlog) + else: + spec.newconst = spec.newconst + CONST_CHANGE * direction + self.compute_potential(r_number, lu_idx) diff --git a/src/disslucc/executors/__init__.py b/src/disslucc/executors/__init__.py index e4d73ae..0ce73a6 100644 --- a/src/disslucc/executors/__init__.py +++ b/src/disslucc/executors/__init__.py @@ -21,8 +21,9 @@ from .base import LuccExecutorBase as LuccExecutorBase from .continuous import LuccContinuousExecutor as LuccContinuousExecutor from .discrete import LuccDiscreteExecutor as LuccDiscreteExecutor + from .saturation import LuccSaturationExecutor as LuccSaturationExecutor -__all__ = ["LuccContinuousExecutor", "LuccDiscreteExecutor", "LuccExecutorBase"] +__all__ = ["LuccContinuousExecutor", "LuccDiscreteExecutor", "LuccExecutorBase", "LuccSaturationExecutor"] def __getattr__(name: str): @@ -35,4 +36,7 @@ def __getattr__(name: str): if name == "LuccDiscreteExecutor": from .discrete import LuccDiscreteExecutor return LuccDiscreteExecutor + if name == "LuccSaturationExecutor": + from .saturation import LuccSaturationExecutor + return LuccSaturationExecutor raise AttributeError(f"module {__name__!r} has no attribute {name!r}") diff --git a/src/disslucc/executors/base.py b/src/disslucc/executors/base.py index 70079c3..3b747e3 100644 --- a/src/disslucc/executors/base.py +++ b/src/disslucc/executors/base.py @@ -44,6 +44,8 @@ from dissmodel.io.convert import vector_to_raster_backend from dissmodel.io.raster import save_geotiff +RASTER_SUFFIXES = (".tif", ".tiff") + class LuccExecutorBase(ModelExecutor): """ @@ -59,19 +61,34 @@ class LuccExecutorBase(ModelExecutor): def validate(self, record: ExperimentRecord) -> None: if not record.source or not record.source.uri: - raise ValueError("record.source.uri is empty -- point it at the input shapefile") + raise ValueError("record.source.uri is empty -- point it at the input GeoTIFF or shapefile") missing = [p for p in self.required_parameters if p not in record.parameters] if missing: raise ValueError(f"Missing parameters: {missing}") def load(self, record: ExperimentRecord) -> RasterBackend: - """Loads the input vector AND ALREADY RASTERIZES -- returns the - RasterBackend ready to use, not the GeoDataFrame. `run()` just - runs the model on top of whatever arrives here.""" + """Returns the RasterBackend ready to use. A GeoTIFF whose bands are + named (a `name` tag per band, as `dissmodel.io.save_geotiff` writes + them) is read as it is -- the land uses, the drivers and any other + band (`mask`, regions, cell order). A vector is loaded AND ALREADY + RASTERIZED, not returned as a GeoDataFrame. `run()` just runs the + model on top of whatever arrives here.""" p = record.parameters land_use_types = p["land_use_types"] driver_names = {k for spec in p["potential_data"] for k in spec.get("betas", {})} + uri = record.source.uri + if uri.lower().endswith(RASTER_SUFFIXES): + (backend, meta), checksum = load_dataset(uri, fmt="raster") + record.source.checksum = checksum + backend.transform = meta["transform"] + backend.crs = meta["crs"] + missing = sorted((set(land_use_types) | driver_names) - set(backend.arrays)) + if missing: + raise ValueError(f"{uri}: bands missing: {missing}") + record.add_log(f"Loaded raster: shape={backend.shape}, {len(backend.arrays)} bands") + return backend + gdf, checksum = load_dataset(record.source.uri, fmt="vector") record.source.checksum = checksum if record.column_map: diff --git a/src/disslucc/executors/saturation.py b/src/disslucc/executors/saturation.py new file mode 100644 index 0000000..4940c6f --- /dev/null +++ b/src/disslucc/executors/saturation.py @@ -0,0 +1,186 @@ +""" +disslucc.executors.saturation +----------------------------- +`LuccSaturationExecutor` -- the continuous model of LuccME-BR (Bezerra et +al. 2022) as an executor: `DemandPreComputedValues` + +`PotentialSpatialLagRegression` + `AllocationClueLikeSaturation`, with +parameters per class and per region, from a model TOML +(`examples/dissmodel-configs/lucc_saturation.toml`). + +Differs from `LuccContinuousExecutor` in three ways, each because of what +this model needs: + +- **raster input.** `record.source.uri` is usually a GeoTIFF whose bands + are named (read by the base's `load()`): the land uses, the drivers and, + optionally, `mask`, `region`, `regionAloc` and the cell order. A + continental cellular space is built once as a raster; rasterizing a + vector at every run would be the wrong way round. +- **regions.** Each `potential_data`/`allocation_data` entry names its + class (`lu`) and, optionally, its `region` (default 1); the executor + groups them into the `[region][class]` lists the components take. +- **intermediate years.** `save_steps` lists steps to write besides the + last, as `_step.tif`, to compare a run against maps of + intermediate years. +""" +from __future__ import annotations + +import pathlib +import time +from typing import Any, ClassVar + +import numpy as np +from dissmodel.core import Environment, Model +from dissmodel.executor import ExperimentRecord +from dissmodel.geo import RasterBackend +from dissmodel.io.raster import save_geotiff + +from ..components.allocation import AllocationClueLikeSaturation +from ..components.demand import DemandPreComputedValues, load_demand_csv +from ..components.potential import PotentialSpatialLagRegression +from ..schemas import SaturationAllocationSpec, SpatialLagRegressionSpec +from .base import LuccExecutorBase + + +def _by_region(entries: list[dict[str, Any]], land_use_types: list[str], build) -> list[list[Any]]: + """[{lu, region?, ...}, ...] -> [[spec of each class, in land_use_types order] per region].""" + regions = sorted({int(e.get("region", 1)) for e in entries}) + if regions != list(range(1, len(regions) + 1)): + raise ValueError(f"regions must be 1..n, got {regions}") + out = [] + for r in regions: + of_region = {e["lu"]: e for e in entries if int(e.get("region", 1)) == r} + missing = [lu for lu in land_use_types if lu not in of_region] + if missing: + raise ValueError(f"region {r}: no entry for {missing}") + out.append([build(of_region[lu]) for lu in land_use_types]) + return out + + +def _potential_spec(d: dict[str, Any]) -> SpatialLagRegressionSpec: + return SpatialLagRegressionSpec( + const=d["const"], ro=d["ro"], betas=dict(d.get("betas", {})), is_log=d.get("is_log", False), + min_reg=d.get("min_reg", 0.0), max_reg=d.get("max_reg", 1.0), + ) + + +def _allocation_spec(d: dict[str, Any]) -> SaturationAllocationSpec: + return SaturationAllocationSpec( + static=d.get("static", -1), + min_value=d.get("min_value", 0.0), max_value=d.get("max_value", 1.0), + min_change=d.get("min_change", 0.0), max_change=d.get("max_change", 1.0), + change_limiar_value=d.get("change_limiar_value", 1.0), + max_change_above_limiar=d.get("max_change_above_limiar", 0.0), + ) + + +class LuccSaturationExecutor(LuccExecutorBase): + """ + Parameters expected in `record.parameters` + --------------------------------------------- + land_use_types : list[str] + complementar_lu : str + demand_csv : str -- path to the demand CSV (columns = land_use_types), one row per step + potential_data : list[dict] -- {"lu", "region"?, "const", "ro", "betas", "is_log"?, "min_reg"?, "max_reg"?} + allocation_data : list[dict] -- {"lu", "region"?, "static", "min_value", "max_value", "min_change", + "max_change", "change_limiar_value", "max_change_above_limiar"} + land_use_no_data, attr_protection, saturation_indicator, order_attr (optional) + cell_area, n_steps, max_difference, max_iteration, initial/min/max_elasticity (optional) + save_steps : list[int] (optional) -- steps written besides the last + + `record.source.uri`: a GeoTIFF with named bands (see module docstring), or + a vector file rasterized by the base (`resolution`). + """ + + name = "lucc_continuous_saturation" + required_parameters: ClassVar[list[str]] = [ + "land_use_types", "demand_csv", "potential_data", "complementar_lu", "allocation_data", + ] + + def run(self, data: RasterBackend, record: ExperimentRecord) -> dict: + p = record.parameters + land_use_types = p["land_use_types"] + backend = data + annual_demand = load_demand_csv(pathlib.Path(p["demand_csv"]).read_text(), land_use_types) + n_steps = int(p.get("n_steps", len(annual_demand))) + cell_area = float(p.get("cell_area", 1.0)) + save_steps = sorted({int(s) for s in p.get("save_steps", [])}) + if "mask" not in backend.arrays: + total = np.sum([np.asarray(backend.get(lu), dtype=np.float64) for lu in land_use_types], axis=0) + backend.set("mask", (total > 0).astype(np.float32)) + + env = Environment(end_time=n_steps - 1) + demand = DemandPreComputedValues(annual_demand=annual_demand, land_use_types=land_use_types) + potential = PotentialSpatialLagRegression( + backend=backend, demand=demand, land_use_types=land_use_types, + potential_data=_by_region(p["potential_data"], land_use_types, _potential_spec), + land_use_no_data=p.get("land_use_no_data"), + ) + allocation = AllocationClueLikeSaturation( + backend=backend, demand=demand, potential=potential, land_use_types=land_use_types, + allocation_data=_by_region(p["allocation_data"], land_use_types, _allocation_spec), + complementar_lu=p["complementar_lu"], cell_area=cell_area, + land_use_no_data=p.get("land_use_no_data"), + attr_protection=p.get("attr_protection"), + saturation_indicator=p.get("saturation_indicator", "saturationLimiar"), + max_difference=float(p.get("max_difference", 1643)), + max_iteration=int(p.get("max_iteration", 1000)), + initial_elasticity=float(p.get("initial_elasticity", 0.1)), + min_elasticity=float(p.get("min_elasticity", 0.001)), + max_elasticity=float(p.get("max_elasticity", 1.5)), + order_attr=p.get("order_attr"), + ) + + snapshots: dict[int, dict[str, np.ndarray]] = {} + + class Snapshot(Model): + def execute(self) -> None: + step = int(self.env.now()) + if step in save_steps: + snapshots[step] = {lu: backend.get(lu).copy() for lu in land_use_types} + + Snapshot() + t0 = time.perf_counter() + env.run() + seconds = time.perf_counter() - t0 + record.add_log(f"Ran {n_steps} steps ({seconds:.1f} s); iterations per step: {allocation.iterations_per_step}") + + mask = backend.get("mask").astype(bool) + areas = {lu: float(backend.get(lu)[mask].sum()) * cell_area for lu in land_use_types} + return { + "backend": backend, + "land_use_types": land_use_types, + "snapshots": snapshots, + "metrics": {f"final_{lu}_area": a for lu, a in areas.items()} + | { + "seconds": seconds, + "iterations_per_step": allocation.iterations_per_step, + "max_error_per_step": allocation.max_error_per_step, + }, + "final_log": f"Final areas: {areas}", + } + + def save(self, result: dict, record: ExperimentRecord) -> ExperimentRecord: + """The base's save() for the last step, plus one GeoTIFF per `save_steps` entry.""" + record = super().save(result, record) + backend, land_use_types = result["backend"], result["land_use_types"] + stem = pathlib.Path(record.output_path) + for step, arrays in sorted(result["snapshots"].items()): + snap = RasterBackend(shape=backend.shape, transform=backend.transform, crs=backend.crs) + for lu, arr in arrays.items(): + snap.set(lu, arr) + snap.set("mask", backend.get("mask")) + uri = str(stem.with_name(f"{stem.stem}_step{step}{stem.suffix or '.tif'}")) + band_spec = [(lu, str(arrays[lu].dtype), -1.0) for lu in land_use_types] + [ + ("mask", str(backend.get("mask").dtype), 0.0) + ] + record.artifacts[f"step{step}"] = save_geotiff( + (snap, {"crs": backend.crs, "transform": backend.transform}), uri, band_spec=band_spec, + ) + record.add_log(f"Saved step {step} to {uri}") + return record + + +if __name__ == "__main__": + # python -m disslucc.executors.saturation run --toml ... --input cellspace.tif --param demand_csv=... + from dissmodel.executor.cli import run_cli + run_cli(LuccSaturationExecutor) diff --git a/src/disslucc/protocols.py b/src/disslucc/protocols.py index f4a4388..9f738e3 100644 --- a/src/disslucc/protocols.py +++ b/src/disslucc/protocols.py @@ -27,3 +27,13 @@ def change_lu_direction(self, lu_index: int) -> int: ... @runtime_checkable class PotentialProtocol(Protocol): def modify(self, r_number: int, lu_idx: int, direction: int) -> None: ... + + +@runtime_checkable +class RegionalPotentialProtocol(PotentialProtocol, Protocol): + """What AllocationClueLikeSaturation needs from its potential: the + number of potential regions (it adapts the elasticities once per + region, as LuccME does) and LuccME's `modifyDriver`, called when the + allocation has not converged after half its iterations.""" + potential_data: list + def modify_driver(self, attr_protection: str, rate: float) -> None: ... diff --git a/src/disslucc/schemas.py b/src/disslucc/schemas.py index b0a5e26..8be9603 100644 --- a/src/disslucc/schemas.py +++ b/src/disslucc/schemas.py @@ -59,3 +59,48 @@ class AllocationSpec: max_value: float = 1.0 min_change: float = 0.0 max_change: float = 1.0 + + +@dataclass +class SpatialLagRegressionSpec: + """ + Parameters of a spatial-lag regression for one land use class. + + Port of one entry of LuccME's `PotentialCSpatialLagRegression.potentialData` + (`isLog`, `const`, `minReg`, `maxReg`, `ro`, `betas`). Used by + PotentialSpatialLagRegression (continuous). + + `newconst` is runtime state, as in RegressionSpec: the constant after + the allocation's adjustments within a step (`modify`). Unlike + RegressionSpec, `const` itself also changes during a run -- LuccME's + `adaptRegressionConstants` writes the adapted value back into it every + step, so the adaptation accumulates (the lab03/lab06 goldens reject the + non-cumulative reading). + """ + const: float + ro: float + betas: dict[str, float] = field(default_factory=dict) + is_log: bool = False + min_reg: float = 0.0 + max_reg: float = 1.0 + newconst: float = field(default=0.0, init=False, repr=False, compare=False) + + +@dataclass +class SaturationAllocationSpec: + """ + Allocation constraints of one land use class in one region, for + AllocationClueLikeSaturation. + + AllocationSpec plus `static` per class (LuccME keeps it in + allocationData) and the two saturation parameters: where the cell's + saturation indicator exceeds `change_limiar_value`, the change in the + demand's direction is halved or capped at `max_change_above_limiar`. + """ + static: int = -1 # -1 = follows demand | 0 = free | 1 = fixed + min_value: float = 0.0 + max_value: float = 1.0 + min_change: float = 0.0 + max_change: float = 1.0 + change_limiar_value: float = 1.0 + max_change_above_limiar: float = 0.0 diff --git a/tests/_lab03_helpers.py b/tests/_lab03_helpers.py new file mode 100644 index 0000000..95fee13 --- /dev/null +++ b/tests/_lab03_helpers.py @@ -0,0 +1,253 @@ +""" +tests/_lab03_helpers.py +----------------------- +Shared by the tests of PotentialSpatialLagRegression and +AllocationClueLikeSaturation against the LuccME goldens of lab03 and lab06 +(benchmark/goldens/, from LambdaGeo/terrame-docker v0.1.1). + +Both labs: DemandPreComputedValues + PotentialCSpatialLagRegression + +AllocationCClueLikeSaturation on csAC (data/input/csAC.zip), 2008–2014, the same +coefficients, demand and allocation data (luccme/tests/functional/lab03.lua, +lab06.lua). lab06 also replaces `ti` from csAC_2009 in 2009 (`updateYears`). +""" + +from __future__ import annotations + +import json +import re +from pathlib import Path + +import geopandas as gpd +import numpy as np +import pandas as pd +from dissmodel.core import Environment, Model +from dissmodel.geo.raster.backend import RasterBackend + +from disslucc import ( + AllocationClueLikeSaturation, + DemandPreComputedValues, + PotentialSpatialLagRegression, + SaturationAllocationSpec, + SpatialLagRegressionSpec, +) +from disslucc.components.potential.spatial_lag import spatial_lag_regression + +ROOT = Path(__file__).resolve().parent.parent +GOLDENS = ROOT / "benchmark" / "goldens" +INPUT = ROOT / "data" / "input" +LAND_USES = ["f", "d", "outros"] +NO_DATA = "outros" +YEARS = list(range(2008, 2015)) +# luccme/tests/functional/lab03.lua and lab06.lua (identical in both) +DEMAND = [ + [137878.1691, 19982.62882, 6489.202049], + [137622.2199, 20238.57805, 6489.202049], + [137366.2707, 20494.52729, 6489.202049], + [137110.3214, 20750.47652, 6489.202049], + [136824.6853, 21036.11265, 6489.202049], + [136539.0492, 21321.74879, 6489.202049], + [136253.4130, 21607.38493, 6489.202049], +] +UPDATES = {"lab03": {}, "lab06": {2009: "csAC_2009.zip"}} + + +def specs() -> list[SpatialLagRegressionSpec]: + return [ + SpatialLagRegressionSpec( + const=0.05266679, + ro=0.9124615, + betas={"uc_us": 0.03789872, "uc_pi": 0.04141921, "ti": 0.04455667}, + ), + SpatialLagRegressionSpec( + const=0.01431553, + ro=0.9019253, + betas={ + "assentamen": 0.0443537, + "uc_us": -0.01454847, + "fertilidad": 0.01701601, + "dist_riobr": -0.00000002262071, + }, + ), + SpatialLagRegressionSpec(const=0, ro=0, betas={}), + ] + + +class Lab: + """csAC on a raster by (row, col), the golden, and the dynamic-variable updates.""" + + def __init__(self, name: str): + self.name = name + self.cells = gpd.read_file(INPUT / "csAC.zip") + self.rows = self.cells["row"].astype(int).values + self.cols = self.cells["col"].astype(int).values + self.shape = (self.rows.max() + 1, self.cols.max() + 1) + self.valid = self.grid(np.ones(len(self.cells))) > 0 + self.golden = pd.read_csv(GOLDENS / name / f"{name}.csv.gz") + self.manifest = json.loads((GOLDENS / name / "manifest.json").read_text()) + self.drivers_used = sorted({d for s in specs() for d in s.betas}) + + def grid(self, values) -> np.ndarray: + a = np.zeros(self.shape) + a[self.rows, self.cols] = np.asarray(values, dtype=np.float64) + return a + + def initial(self, column: str) -> np.ndarray: + return self.grid(self.cells[column].astype(float).values) + + def drivers(self, year: int, pair_by: str = "position") -> dict[str, np.ndarray]: + """Drivers in effect in `year`. LuccME's updateYears copies the columns + of the _ file into the cells *by position* (forEachCellPair), + although csAC_2009 is not in the same order as csAC; `pair_by="id"` + pairs by object_id0 instead (used to show the difference matters).""" + cells = self.cells.copy() + for upd_year, file in sorted(UPDATES[self.name].items()): + if year >= upd_year: + new = gpd.read_file(INPUT / file) + if pair_by == "id": + new = cells[["object_id0"]].merge( + new.drop(columns="geometry"), on="object_id0", how="left" + ) + for col in new.columns: + if col in cells.columns and col not in ("geometry", "object_id0"): + cells[col] = new[col].values + return {d: self.grid(cells[d].astype(float).values) for d in self.drivers_used} + + def golden_year(self, year: int) -> pd.DataFrame: + return self.golden[self.golden["year"] == year].set_index(["row", "col"]) + + def past(self, year: int) -> dict[str, np.ndarray]: + """Land use at the start of `year`: initial attributes, then the golden's previous year.""" + if year == YEARS[0]: + return {lu: self.initial(lu) for lu in LAND_USES} + prev = self.golden_year(year - 1).reindex(list(zip(self.rows, self.cols))) + return {lu: self.grid(prev[f"{lu}_out"].values) for lu in LAND_USES} + + +def adapted_consts(cumulative: bool = True) -> dict[int, list[float]]: + """const in effect each year. LuccME writes the adapted constant back into + `const`, so the adaptation accumulates; `cumulative=False` is the other + reading of the code, kept to show the golden tells them apart.""" + base = [s.const for s in specs()] + out, cur = {}, list(base) + for t, year in enumerate(YEARS): + if t > 0: + start = cur if cumulative else base + cur = [ + start[i] + 0.01 * (DEMAND[t][i] - DEMAND[t - 1][i]) / DEMAND[t - 1][i] + for i in range(len(LAND_USES)) + ] + out[year] = list(cur) + return out + + +def max_errors(lab: Lab, cumulative: bool = True, pair_by: str = "position") -> dict[tuple[int, str], float]: + consts = adapted_consts(cumulative) + errors = {} + for year in YEARS: + past = lab.past(year) + drivers = lab.drivers(year, pair_by) + ref = lab.golden_year(year) + r = ref.index.get_level_values("row").values + c = ref.index.get_level_values("col").values + for i, (lu, spec) in enumerate(zip(LAND_USES, specs(), strict=True)): + spec.const = consts[year][i] + _, pot = spatial_lag_regression(spec, past[lu], drivers, lab.valid, past[NO_DATA]) + errors[(year, lu)] = float(np.abs(pot[r, c] - ref[f"{lu}_pot"].values).max()) + return errors + + +CELL_AREA = 25 # lab03.lua / lab06.lua: cellArea = 25 +# lab03.lua / lab06.lua allocationData, region 1 +ALLOCATION = [ + SaturationAllocationSpec(static=-1, min_value=0, max_value=1, min_change=0, max_change=1), # f + SaturationAllocationSpec(static=-1, min_value=0, max_value=1, min_change=0, max_change=1), # d + SaturationAllocationSpec(static=1, min_value=0, max_value=1, min_change=0, max_change=1), # outros +] + + +def terrame_log(name: str) -> dict[int, tuple[int, float]]: + """year → (iterations, maximum error) from the golden's TerraME log.""" + text = (GOLDENS / name / "terrame.log").read_text() + pattern = ( + r"Demand allocated correctly in (\d+)\s+Number of iterations: (\d+)\s+Maximum error: ([\d.eE+-]+)" + ) + return {int(y): (int(n), float(e)) for y, n, e in re.findall(pattern, text)} + + +def run(lab: Lab, order: str = "file"): + """Run the model from 2008 to 2014; return per-year snapshots and the allocation.""" + backend = RasterBackend(shape=lab.shape) + for lu in LAND_USES: + backend.set(lu, lab.initial(lu)) + for name, arr in lab.drivers(YEARS[0]).items(): + backend.set(name, arr) + backend.set("mask", lab.valid.astype(np.float32)) + # the order in which TerraME visits the cells: the layer's (file) order + visit = ( + np.arange(len(lab.cells)) if order == "file" else np.random.default_rng(0).permutation(len(lab.cells)) + ) + backend.set("order", lab.grid(visit)) + + class Updates(Model): + """LuccME's updateYears: this year's drivers, before demand and potential.""" + + def execute(self): + year = YEARS[0] + int(self.env.now()) + for name, arr in lab.drivers(year).items(): + backend.set(name, arr) + + snaps: dict[int, dict[str, np.ndarray]] = {} + + class Record(Model): + def execute(self): + year = YEARS[0] + int(self.env.now()) + snaps[year] = { + **{f"{lu}_out": backend.get(lu).copy() for lu in LAND_USES}, + **{f"{lu}_pot": backend.get(f"{lu}_pot").copy() for lu in LAND_USES}, + } + + env = Environment(end_time=len(YEARS) - 1) + Updates() + demand = DemandPreComputedValues(annual_demand=DEMAND, land_use_types=LAND_USES) + potential = PotentialSpatialLagRegression( + backend=backend, + potential_data=[specs()], + demand=demand, + land_use_types=LAND_USES, + land_use_no_data=NO_DATA, + ) + allocation = AllocationClueLikeSaturation( + backend=backend, + demand=demand, + potential=potential, + land_use_types=LAND_USES, + allocation_data=[ALLOCATION], + complementar_lu="f", + cell_area=CELL_AREA, + land_use_no_data=NO_DATA, + attr_protection="uc_pi", + max_difference=1643, + max_iteration=1000, + initial_elasticity=0.1, + min_elasticity=0.001, + max_elasticity=1.5, + order_attr="order", + ) + Record() + env.run() + return snaps, allocation + + +def worst_per_column(lab: Lab, snaps) -> dict[str, tuple[int, float]]: + worst: dict[str, tuple[int, float]] = {} + for year in YEARS: + ref = lab.golden_year(year) + r = ref.index.get_level_values("row").values + c = ref.index.get_level_values("col").values + for col, arr in snaps[year].items(): + err = float(np.abs(arr[r, c] - ref[col].values).max()) + if col not in worst or err > worst[col][1]: + worst[col] = (year, err) + return worst + + diff --git a/tests/lua/AllocationCClueLikeSaturation.lua b/tests/lua/AllocationCClueLikeSaturation.lua new file mode 100644 index 0000000..7132db8 --- /dev/null +++ b/tests/lua/AllocationCClueLikeSaturation.lua @@ -0,0 +1,858 @@ +--- Modification of the AllocationClueLike component. In this case the speed of change in each cell use a spatio-temporal variable, dynamically updated +-- every year, that indicates if the cell is in a more consolidated or in a frontier area. The saturation threshold considers a 10x10 neighbourhood, +-- discounting protected areas. +-- @arg component.maxDifference Maximum allocation error allowed for each land use in area. +-- @arg component.maxIteration Maximum number of iterations at each time step of the model. +-- @arg component.initialElasticity Initial elasticity value (iterationFactor). +-- @arg component.minElasticity Minimum elasticity value which controls the allocation interaction factor. +-- @arg component.maxElasticity Maximum elasticity value which controls the allocation interaction factor. +-- @arg component.complementarLU The land use which will be recomputed in the end to sum exactly 100%. +-- @arg component.saturationIndicator Name of a attribute which will be dynamically updated (and can be saved for calibration purposes). +-- @arg component.attrProtection Database attribute indicating the percentage of protected areas to be excluded from the saturation level computation. +-- @arg component.allocationData A table with two allocation parameters for each land use. +-- @arg component.allocationData.static Indicates if the variable can increase or decrease in each cell, or only change in the direction of the demand. +-- @arg component.allocationData.minValue Minimum value allowed for the percentage of a given land use in a cell (as a result of new changes - the original +-- value can be out of the limits). +-- @arg component.allocationData.maxValue Maximum value allowed for the percentage of a given land use in a cell (as a result of new changes - the original +-- value can be out of the limits). +-- @arg component.allocationData.minChange Minimum change in a given land use in a cell in a time step until (saturation) threshold. +-- @arg component.allocationData.maxChange Maximum change in a given land use allowed in a cell in a time step until (saturation) threshold. +-- @arg component.allocationData.changeLimiarValue Threshold (or limier) refers to a given amount of the land use in each cell. After this limier, the speed +-- of change of a given land use in the cell is modified. +-- @arg component.allocationData.maxChangeAboveLimiar Maximum change in a given land use allowed in a cell in a time step after (saturation) threshold. +-- @arg component.run Handles with the rules of the component execution. +-- @arg component.verify Handles with the verify method of a AllocationCClueLikeSaturation component. +-- @arg component.updateAllocationParameters Handles with the allocation parameters update. +-- @arg component.initElasticity Handles with the elasticity initialize considering a single +-- elasticity for each land use (all cells). +-- @arg component.computeChange Compute the Allocation change based on the potential of the cell. +-- @arg component.compareAllocationToDemand Compares the demand to the amount of allocated land +-- use/cover, then adapts elasticity. +-- @arg component.correctCellChange Corrects total land use/cover types to 100 percent. +-- @arg component.countAllocatedLandUseArea Calculates total area allocated by the regression +-- equations for each land use/cover type. +-- @arg component.printAllocatedArea Calculates and prints the allocated by the regression +-- equations for each land use/cover type. +-- @usage --DONTRUN +--A1 = AllocationCClueLikeSaturation +--{ +-- maxDifference = 1643, +-- maxIteration = 1000, +-- initialElasticity = 0.1, +-- minElasticity = 0.001, +-- maxElasticity = 1.5, +-- complementarLU = "floresta", +-- saturationIndicator = "saturationLimiar", +-- attrProtection = "uc_pi", +-- allocationData = +-- { +-- -- Region 1 +-- { +-- {static = -1, minValue = 0, maxValue = 1, minChange = 0, maxChange = 1, changeLimiarValue = 1, maxChangeAboveLimiar = 0}, -- floresta +-- {static = -1, minValue = 0, maxValue = 1, minChange = 0, maxChange = 1, changeLimiarValue = 1, maxChangeAboveLimiar = 0}, -- desmatamento +-- {static = 1, minValue = 0, maxValue = 1, minChange = 0, maxChange = 1, changeLimiarValue = 1, maxChangeAboveLimiar = 0}, -- outros +-- } +-- } +--} +function AllocationCClueLikeSaturation(component) + -- Handles with the rules of the component execution. + -- @arg event A representation of a time instant when the simulation engine must run. + -- @arg luccMEModel A LuccME model. + -- @usage --DONTRUN + -- component.run(event, model) + component.run = function(self,event,luccMEModel) + -- Synchronize cellular space in the first year + local luTypes = luccMEModel.landUseTypes + local cs = luccMEModel.cs + + -- saving the first year values + if (event:getTime() == luccMEModel.startTime) then + for k, cell in pairs (cs.cells) do + for luind, lu in pairs (luTypes) do + cell[lu.."_backup"] = cell[lu] + end + end + + -- restoring the data from a saved year + elseif (belong(event:getTime() - 1,luccMEModel.save.saveYears) and ((event:getTime() - 1) ~= luccMEModel.startTime)) then + for k, cell in pairs (cs.cells) do + for luind, lu in pairs (luTypes) do + cell[lu] = cell[lu.."_backupYear"] + end + end + end + + -- Initialize the demandDirection and elasticity(internal component variables) + self:initElasticity(luccMEModel, self.initialElasticity) + self:updateAllocationParameters(event, luccMEModel) + + -- Define iteration loop variables + local nIter = 0 + local allocation_ok = false + local maxAdjust = self.maxDifference + local maxdiff = self.maxDifference * 1000 -- Used to have a large number of iterations + local flagFlex = false + local regionsNumber = #self.allocationData + + -- Loop until maxdiff is achieved + repeat + for rNumber = 1, regionsNumber, 1 do + -- compute tentative allocation + if (event:getTime() ~= luccMEModel.startTime) then + self:computeChange(luccMEModel, rNumber) + self:correctCellChange(luccMEModel, rNumber) + end + end + + if luccMEModel.useLog == true then + self:printAllocatedArea(event, luccMEModel, nIter) + end + + -- verify if allocation reaches demand + maxdiff = self:compareAllocationToDemand(event, luccMEModel) + + if (maxdiff <= maxAdjust) then + allocation_ok = true + + if (luccMEModel.useLog == true) then + print("\nDemand allocated correctly in "..event:getTime().."\tNumber of iterations: "..nIter.."\tMaximum error: "..maxdiff) + end + else + nIter = nIter + 1 + + if (nIter > self.maxIteration * 0.50) and (flagFlex == false) then + maxAdjust = maxAdjust * 2 + flagFlex = true + luccMEModel.potential:modifyDriver(self.complementarLU, self.attrProtection, 0.5, event, luccMEModel) + end + end + until ((nIter >= self.maxIteration) or (allocation_ok == true)) + + if (nIter == self.maxIteration) then + print("\nDemand not allocated correctly in this time step: "..nIter.."\n") + os.exit() + end + + forEachCell(cs, function(cell) + local total = 0 + local biggerLuValue = 0 + local biggerLu = "" + + for i, lu in pairs (luTypes) do + if (lu ~= self.complementarLU) then + total = total + cell[lu] + if (lu ~= luccMEModel.landUseNoData) then + if (cell[lu] > biggerLuValue) then + biggerLuValue = cell[lu] + biggerLu = lu + end + end + end + end + if (self.complementarLU ~= nil) then + cell[self.complementarLU] = 1 - total + if (cell[self.complementarLU] < 0) then + cell[biggerLu] = cell[biggerLu] + cell[self.complementarLU] + cell[self.complementarLU] = 0 + end + end + end + ) + + -- calculating changes and outputs + for i, lu in pairs (luTypes) do + local out = lu.."_out" + local diff = lu.."_chpast" + local previous = lu.."_chtot" + + for k, cell in pairs (cs.cells) do + if (event:getTime() == luccMEModel.startTime) then + cell[previous] = 0 + else + cell[previous] = cell[lu] - cell[lu.."_backup"] + end + cell[out] = cell[lu] + cell[diff] = cell[lu] - cell.past[lu] + end + end + + -- doing save year backup + if (belong(event:getTime(),luccMEModel.save.saveYears) and (event:getTime() ~= luccMEModel.startTime)) then + for k, cell in pairs (cs.cells) do + for luind, lu in pairs (luTypes) do + cell[lu.."_backupYear"] = cell[lu] + cell[lu] = cell[lu.."_backup"] + end + end + end + + cs:synchronize() + end + + -- Handles with the parameters verification. + -- @arg event A representation of a time instant when the simulation engine must run. + -- @arg luccMEModel A LuccME model. + -- @usage --DONTRUN + -- component.verify(event, self) + component.verify = function(self, event, luccMEModel) + print("Verifying Allocation parameters") + + -- create a region if does not have one + local cs = luccMEModel.cs + + if (self.regionAttr == nil) then + self.regionAttr = "regionAloc" + end + + forEachCell(cs, function(cell) + cell["alternate_model"] = 0 + + if (cell[self.regionAttr] == nil) then + cell["regionAloc"] = 1 + else + cell["regionAloc"] = cell[self.regionAttr] + end + end + ) + + -- check maxDifference + if (self.maxDifference == nil) then + error("maxDifference variable is missing", 2) + end + + -- check maxIteration + if (self.maxIteration == nil) then + error("maxIteration variable is missing", 2) + end + + -- check initialElasticity + if (self.initialElasticity == nil) then + error("initialElasticity variable is missing", 2) + end + + -- check minElasticity + if (self.minElasticity == nil) then + error("minElasticity variable is missing", 2) + end + + -- check maxElasticity + if (self.maxElasticity == nil) then + error("maxElasticity variable is missing", 2) + end + + -- check complementarLU + if (self.complementarLU == nil) then + error("complementarLU variable is missing", 2) + end + + -- check saturationIndicator + if (self.saturationIndicator == nil) then + error("saturationIndicator variable is missing", 2) + end + + -- check attrProtection + if (self.attrProtection == nil) then + error("attrProtection variable is missing", 2) + end + + -- check allocationData + if (self.allocationData == nil) then + error("allocationData is missing", 2) + end + + local regionsNumber = #self.allocationData + local lutNumber = #luccMEModel.landUseTypes + + -- check number of Regions + if (regionsNumber == nil or regionsNumber == 0) then + error("The model must have at least One region") + else + for i = 1, regionsNumber, 1 do + local allocationNumber = #self.allocationData[i] + + -- check the number of allocations + if (allocationNumber ~= lutNumber) then + error("Invalid number of regressions on Region number "..i.." . Regressions: "..allocationNumber.." LandUseTypes: "..lutNumber) + end + + -- check data into allocationData + for j = 1, allocationNumber, 1 do + -- check static variable + if(self.allocationData[i][j].static == nil) then + error("static variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + + -- check minValue variable + if (self.allocationData[i][j].minValue == nil) then + error("minValue variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + + -- check maxValue variable + if (self.allocationData[i][j].maxValue == nil) then + error("maxValue variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + + -- check minChange variable + if (self.allocationData[i][j].minChange == nil) then + error("minChange variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + + -- check maxChange variable + if (self.allocationData[i][j].maxChange == nil) then + error("maxChange variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + + -- check changeLimiarValue variable + if (self.allocationData[i][j].changeLimiarValue == nil) then + error("changeLimiarValue variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + + -- check maxChangeAboveLimiar variable + if (self.allocationData[i][j].maxChangeAboveLimiar == nil) then + error("maxChangeAboveLimiar variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + end -- for j + end -- for i + end -- else + + -- check complementarLU within database + if (self.complementarLU ~= nil) then + local foundCLU = 0 + + for k = 1, lutNumber, 1 do + if (self.complementarLU == luccMEModel.landUseTypes[k]) then + foundCLU = 1 + break + end + end + + if (foundCLU == 0) then + error("complementarLU: "..self.complementarLU.." is not a landUseType defined on main file.", 2) + end + end + + -- check attrProtection within database + if (self.attrProtection ~= nil) then + if (luccMEModel.cs.cells[1][self.attrProtection] == nil) then + error("attrProtection: "..self.attrProtection.." not found within database", 2) + end + end + + luccMEModel.cs:createNeighborhood{ name = "11x11", + strategy = "mxn", + } + end + + -- Update the allocation parameters based on the saturation of the region. + -- @arg event A representation of a time instant when the simulation engine must run. + -- @arg luccMEModel A LuccME model. + -- @usage --DONTRUN + -- component.updateAllocationParameters(event, luccMEModel) + component.updateAllocationParameters = function(self, event, luccMEModel) + local cs = luccMEModel.cs + local currentTime = event:getTime() + + forEachCell(cs, function(cell) + local prot_t = 0 + local original = 1 + local perc_def_original = 0 + + if (self.attrProtection ~= nil) then + prot_t = cell[self.attrProtection] + end + + if (luccMEModel.landUseNoData ~= nil) then + original = 1 - cell[luccMEModel.landUseNoData] + end + + local available_forest = original - prot_t + + if (available_forest > 0) then + local total_perc = 0 + local count = 0 + + perc_def_original = (1 - cell[self.complementarLU]) / available_forest + + if (perc_def_original > 1) then + perc_def_original = 1 + end + + forEachNeighbor(cell, "11x11", function(neigh, _, cell) + local prot_t = 0 + local original = 1 + + if (self.attrProtection ~= nil) then + prot_t = neigh[self.attrProtection] + end + + if (luccMEModel.landUseNoData ~= nil) then + original = 1 - neigh[luccMEModel.landUseNoData] + end + + local neigh_available_forest = original - prot_t + + if (neigh_available_forest > 0) then + local neigh_perc = (1 - neigh[self.complementarLU]) / neigh_available_forest + if (neigh_perc > 1) then + neigh_perc = 1 + end + + total_perc = total_perc + neigh_perc + count = count + 1 + end + end + ) + + if (count > 0) then + perc_def_original = total_perc / count + end + else + perc_def_original = 1 + end + + cell[self.saturationIndicator] = perc_def_original + end + ) + end + + -- Handles with the elasticity initialize considering a single elasticity for each land use (all cells). + -- Similar to the coarse scale old clue. + -- @arg luccMEModel A LuccME model. + -- @arg value The elasticity value. + -- @usage --DONTRUN + -- component.initElasticity(luccMEModel, self.initialElasticity) + component.initElasticity = function(self, luccMEModel, value) + -- Init elasticity. In this version of the component, a single elasticity for each land use(all cells). + -- Similar to the coarse scale old clue + local luTypes = luccMEModel.landUseTypes + self.elasticity = {} + + for k = 1, #luTypes, 1 do + self.elasticity[k] = value + end + end + + -- Compute the allocation change. + -- @arg luccMEModel A LuccME model. + -- @usage --DONTRUN + -- component.computeChange(luccMEModel) + component.computeChange = function(self, luccMEModel, rNumber) + local cs = luccMEModel.cs + local luTypes = luccMEModel.landUseTypes + local activeRegionNumber = 0 + + for i, luAllocData in pairs (self.allocationData[rNumber]) do + local lu = luTypes[i] + local attr_pot = lu.."_pot" + local luDirect = luccMEModel.demand:getCurrentLuDirection(i) + + for k, cell in pairs (cs.cells) do + if (cell.regionAloc == rNumber) then + local pot = cell[attr_pot] + local luStatic = luAllocData.static + local change = pot * self.elasticity[i] + + activeRegionNumber = rNumber + + if (math.abs(change) < luAllocData.minChange) then + pot = 0 + change = 0 + end + + if (math.abs(change) >= luAllocData.maxChange) then + if (pot ~= 0) then + change = luAllocData.maxChange * (pot / math.abs(pot)) + end + end + + if (((pot >= 0) and (luDirect == 1) and (luStatic < 1) and (cell[self.saturationIndicator] > luAllocData.changeLimiarValue))) then + if (change >= luAllocData.maxChangeAboveLimiar) then + if ((change / 2) < luAllocData.maxChangeAboveLimiar) then + change = change / 2 + else + change = luAllocData.maxChangeAboveLimiar + end + end + end + + if ((pot <= 0) and (luDirect == -1) and (luStatic < 1) and (cell[self.saturationIndicator] > luAllocData.changeLimiarValue)) then + if (math.abs(change) >= luAllocData.maxChangeAboveLimiar) then + if (math.abs(change / 2) < luAllocData.maxChangeAboveLimiar) then + change = change / 2 + else + change = (-1) * luAllocData.maxChangeAboveLimiar + end + end + end + + + if (luStatic == 1) then --do not change + cell[lu] = cell.past[lu] + elseif (luStatic == 0) then -- change independent of demand direction (ANAP) + cell[lu] = cell.past[lu] + change + elseif (((pot >= 0) and (luDirect == 1) and (luStatic == -1)) or ((pot <= 0) and (luDirect == -1) and (luStatic == -1))) then + cell[lu] = cell.past[lu] + change -- change according to demand direction + else + cell[lu] = cell.past[lu] + end + + if (cell[lu] < 0) then + cell[lu] = 0 + end + + if (cell[lu] < luAllocData.minValue) then + if (cell.past[lu] >= luAllocData.minValue) then + cell[lu] = luAllocData.minValue + else + cell[lu] = cell.past[lu] + end + end + + if (cell[lu] > luAllocData.maxValue) then + if (cell.past[lu] <= luAllocData.maxValue) then + cell[lu] = luAllocData.maxValue + else + cell[lu] = cell.past[lu] + end + end + end -- if cell.region + end -- for k + + if (activeRegionNumber == 0) then + error("Region "..rNumber.." is not set into database.") + end + end -- for i, lu + end -- computeChange + + -- Compares the demand to the amount of allocated land use/cover, then adapts elasticity. + -- @arg event A representation of a time instant when the simulation engine must run. + -- @arg luccMEModel A LuccME model. + -- @usage --DONTRUN + -- component.compareAllocationToDemand (event, luccMEModel) + component.compareAllocationToDemand = function(self, event, luccMEModel, rNumber) + -- Compares the demand to the amount of allocated land use/cover, then adapts elasticity + local cs = luccMEModel.cs + local luTypes = luccMEModel.landUseTypes + local areas = self:countAllocatedLandUseArea(cs, luTypes) + local max = 0 + local tot = 0 + local nRegression = 0 + local regionsNumber = #self.allocationData + + for j = 1, #luccMEModel.potential.potentialData, 1 do + for i, lu in pairs (luTypes) do + local luDirect = luccMEModel.demand:getCurrentLuDirection(i) + local currentDemand = luccMEModel.demand:getCurrentLuDemand(i) + + if (luDirect == 0 and event:getTime() == luccMEModel.startTime) then + if (currentDemand >= areas[i]) then + luDirect = 1 + else + luDirect = -1 + end + end + + if (luDirect == 1) then + self.elasticity[i] = self.elasticity[i] * (currentDemand / areas[i]) + else + self.elasticity[i] = self.elasticity[i] * (areas[i] / currentDemand) + end + + if (self.elasticity[i] > self.maxElasticity) then + self.elasticity[i] = self.maxElasticity + luccMEModel.potential:modify(luccMEModel, j, i, luDirect, event) + end + + if (self.elasticity[i] < self.minElasticity) then + if (self.allocationData[1][i].static < 0) then + self.elasticity[i] = self.minElasticity + luccMEModel.potential:modify(luccMEModel, j, i, luDirect * (-1), event) -- Original clue does not modify in this case, but AMAZALERT results are like this + else + luccMEModel.demand:changeLuDirection(i) + end + end + + if (luccMEModel.useLog == true) then + if (j > nRegression) then + print("Region "..j) + nRegression = j + end + + if (luccMEModel.potential.potentialData[j][i].newminReg ~= nil) then + print(lu, "elas: ", self.elasticity[i], "dir: ", luDirect, "const: ", luccMEModel.potential.potentialData[j][i].const, "->", luccMEModel.potential.potentialData[j][i].newconst, luccMEModel.potential.potentialData[j][i].newminReg, luccMEModel.potential.potentialData[j][i].newmaxReg) + elseif (luccMEModel.potential.potentialData[j][i].const ~= nil) then + print(lu, "elas: ", self.elasticity[i], "dir: ", luDirect, "const: ", luccMEModel.potential.potentialData[j][i].const, "->", luccMEModel.potential.potentialData[j][i].newconst) + else + print(lu, "elas: ", self.elasticity[i], "dir: ", luDirect) + end + end + + local diff = math.abs((areas[i] - currentDemand)) + + if (diff > max) then + max = diff + end + + tot = tot + math.abs(areas[i] - currentDemand) + end -- for i + end -- for j + + return max + end + + -- Corrects total land use/cover types to 100 percent. + -- @arg luccMEModel A LuccME model. + -- @usage --DONTRUN + -- component.correctCellChange(luccMEModel) + component.correctCellChange = function(self, luccMEModel, rNumber) + -- corrects total land use/cover types to 100 percent + local cs = luccMEModel.cs + local luTypes = luccMEModel.landUseTypes + local NCOV = #luTypes + local BACKP = 0.5 + + for k, cell in pairs (cs.cells) do + if (cell.regionAloc == rNumber) then + local nostatic = 0 + local decr = 0 + local incr = 0 + local totcov = 0 + local totchange = 0 + local totstatic = 0 + local amin = cell[luTypes[1]] - cell.past[luTypes[1]] + local amax = amin + local max = math.abs(amax) + local l = 0 + + -- checks if total land use/covers from 100 percent + for i, lu in pairs (luTypes) do + totcov = totcov + cell[lu] + end + + if (math.abs(totcov - 1) > 0.005) then + for i, lu in pairs (luTypes) do + local dif = cell[lu] - cell.past[lu] + totchange = totchange + math.abs(dif) + + if (math.abs(dif) > max) then + max = math.abs(dif) + end + + if (dif > amax) then + amax = dif + end + + if (dif < amin) then + amin = dif + end + end + + -- adapts land use/cover types if all of them change into the same direction + if (totchange > 0) then + if ((BACKP * totchange) > (max * 0.5)) then + BACKP =(max / (2 * totchange)) + end + end + + for i, luAllocData in pairs (self.allocationData[rNumber]) do + local lu = luTypes[i] + local luStatic = luAllocData.static + + if ((cell[lu] <= luAllocData.minValue) or cell[lu] >= luAllocData.maxValue) then + luStatic = 1 + end + + if (luStatic < 1) then + local dif = cell[lu] - cell.past[lu] + + if (dif >(-1 * BACKP * totchange)) then + incr = incr + 1 + end + + if (dif <(BACKP * totchange)) then + decr = decr + 1 + end + else + nostatic = nostatic + 1 + end + end + + if ((decr == (NCOV - nostatic)) or (incr == (NCOV - nostatic))) then + for i, luAllocData in pairs (self.allocationData[rNumber]) do + local lu = luTypes[i] + local luDirect = luccMEModel.demand:getCurrentLuDirection(i) + + local luStatic = luAllocData.static + + if ((cell[lu] <= luAllocData.minValue) or cell[lu] >= luAllocData.maxValue) then + luStatic = 1 + end + + if (luStatic < 1) then + if (incr == (NCOV - nostatic)) then + cell[lu] = cell[lu] - (amin +(BACKP * totchange)) + end + + if (decr == (NCOV - nostatic)) then + cell[lu] = cell[lu] + ((BACKP * totchange) - amax) + end + + if (cell[lu] < 0) then + cell[lu] = 0 + end + end + + if ((luDirect == 1) and (cell[lu] < cell.past[lu]) and (luAllocData.static == -1)) then + cell[lu] = cell.past[lu] + end + + if ((luDirect == -1) and (cell[lu] > cell.past[lu]) and (luAllocData.static == -1)) then + cell[lu] = cell.past[lu] + end + end + end + + -- perform the corrections + repeat + l = l + 1 + totcov = 0 + totchange = 0 + totstatic = 0 + + for i, luAllocData in pairs (self.allocationData[rNumber]) do + local lu = luTypes[i] + local luStatic = luAllocData.static + + if ((cell[lu] <= luAllocData.minValue) or cell[lu] >= luAllocData.maxValue) then + luStatic = 1 + end + + if (luStatic < 1) then + totcov = totcov + cell[lu] + else + totstatic = totstatic + cell[lu] + end + + totchange = totchange + math.abs(cell[lu] - cell.past[lu]) + end + + if (math.abs(totcov - (1 - totstatic)) > 0.005) then + if (totchange == 0) then + for i, luAllocData in pairs (self.allocationData[rNumber]) do + local lu = luTypes[i] + local luStatic = luAllocData.static + + if ((cell[lu] <= luAllocData.minValue) or cell[lu] >= luAllocData.maxValue) then + luStatic = 1 + end + + if (luStatic < 1) then + cell[lu] = cell[lu] * ((1 - totstatic) / totcov) + end + end + else + for i, luAllocData in pairs (self.allocationData[rNumber]) do + local lu = luTypes[i] + local aux = cell[lu] + + cell[lu] = cell[lu] -(math.abs(cell[lu] - cell.past[lu]) * ((totcov - (1 - totstatic)) / totchange)) + + if (cell[lu] < 0) then + cell[lu] = 0 + end + end + end + end + until ((math.abs(totcov -(1-totstatic)) <= 0.005) or (l >= 25)) + + if (l == 25) then + totcov = 0 + totstatic = 0 + + for i, luAllocData in pairs (self.allocationData[rNumber]) do + local lu = luTypes[i] + local luStatic = luAllocData.static + + if ((cell[lu] <= luAllocData.minValue) or (cell[lu] >= luAllocData.maxValue)) then + luStatic =1 + end + + if (luStatic < 1) then + totcov = totcov + cell[lu] + else + totstatic = totstatic + cell[lu] + end + end + + for i, luAllocData in pairs (self.allocationData[rNumber]) do + local lu = luTypes[i] + local luStatic = luAllocData.static + + if ((cell[lu] <= luAllocData.minValue) or (cell[lu] >= luAllocData.maxValue)) then + luStatic = 1 + end + + if (luStatic < 1) then + cell[lu] = cell[lu] * ((1 - totstatic) / totcov) + end + end + end -- if l + end -- totcov + end -- if cell.region + end -- for each cell + end + + -- Calculates total area allocated by the regression equations for each land use/cover type. + -- @arg cs A multivalued set of Cells (Cell Space). + -- @arg luTypes A set of land use types. + -- @usage --DONTRUN + -- component.countAllocatedLandUseArea(cs, luTypes) + component.countAllocatedLandUseArea = function(self, cs, luTypes) + -- Calculates total area allocated by the regression equations for each land use/cover type + local areas = {} + local cellarea = cs.cellArea + + for i, lu in pairs (luTypes) do + local area = 0 + + for k, cell in pairs (cs.cells) do + local temp = cell[lu] + + if (temp > 0) then + area = area + temp + end + end + + areas[i] = area * cellarea + end + + return areas + end + + -- Calculates and prints the allocated by the regression equations for each land use/cover type. + -- @arg event A representation of a time instant when the simulation engine must run. + -- @arg luccMEModel A LuccME model. + -- @arg nIter An iterator number. + -- @usage --DONTRUN + -- component.printAllocatedArea(event, luccMEModel, nIter) + component.printAllocatedArea = function(self, event, luccMEModel, nIter) + -- Calculates and prints the allocated by the regression equations for each land use/cover type + local cs = luccMEModel.cs + local luTypes = luccMEModel.landUseTypes + local areas = {} + + areas = self:countAllocatedLandUseArea(cs, luTypes) + + print("\nYear: "..event:getTime().."\tStep: "..nIter) + + for i, lu in pairs (luTypes) do + local currentDemand = luccMEModel.demand:getCurrentLuDemand(i) + print(lu.."\tarea: "..math.floor(areas[i]).."\tDifference: "..math.floor(areas[i] - currentDemand)) + end + + io.flush() + end + + collectgarbage("collect") + return component +end -- close Allocation Component diff --git a/tests/lua/PotentialCSpatialLagRegression.lua b/tests/lua/PotentialCSpatialLagRegression.lua new file mode 100644 index 0000000..6db3ca3 --- /dev/null +++ b/tests/lua/PotentialCSpatialLagRegression.lua @@ -0,0 +1,411 @@ +--- Similar to the LinearRegression approach, but relies spatial regression techniques to estimate the +-- regression cover (considers the spatial dependence of the land use). +-- @arg component A Spatial Lag Regression component. +-- @arg component.potentialData A table with the regression parameters for each attribute. +-- @arg component.potentialData.isLog Inform whether the model is part of a coupling model. +-- @arg component.potentialData.const A linear regression constant. +-- @arg component.potentialData.minReg A coefficient to minimize the regression value. +-- @arg component.potentialData.maxReg A coefficient to potentiality the regression value. +-- @arg component.potentialData.ro Auto regressive coefficient. +-- @arg component.potentialData.betas A linear regression betas for land use drivers +-- and the index of landUseDrivers to be used by the regression (attributes). +-- @arg component.landUseDrivers The land use drivers fields in database. +-- @arg component.run Handles with the execution method of a PotentialCSpatialLagRegression component. +-- @arg component.verify Handles with the verify method of a PotentialCSpatialLagRegression component. +-- @arg component.modify Handles with the modify method of a PotentialCSpatialLagRegression component. +-- @arg component.adaptRegressionConstants Handles with the constants regression method of a +-- PotentialCSpatialLagRegression component. +-- @arg component.modifyDriver Modify potencial for an protected area. +-- @arg component.computePotential Handles with the modify method of a PotentialCSpatialLagRegression component. +-- @return The modified component. +-- @usage --DONTRUN +--P1 = PotentialCSpatialLagRegression +--{ +-- potentialData = +-- { +-- -- Region 1 +-- { +-- -- floresta +-- { +-- isLog = false, +-- const = 0.05266679, +-- minReg = 0, +-- maxReg = 1, +-- ro = 0.9124615, +-- +-- betas = +-- { +-- uc_us = 0.03789872, +-- uc_pi = 0.04141921, +-- ti = 0.04455667 +-- } +-- }, +-- +-- -- desmatamento +-- { +-- isLog = false, +-- const = 0.01431553, +-- minReg = 0, +-- maxReg = 1, +-- ro = 0.9019253, +-- +-- betas = +-- { +-- assentamentos = 0.0443537, +-- uc_us = -0.01454847, +-- dist_riobranco = -0.00000002262071, +-- fertilidadealtaoumedia = 0.01701601 +-- } +-- }, +-- +-- -- outros +-- { +-- isLog = false, +-- const = 0, +-- minReg = 0, +-- maxReg = 1, +-- ro = 0, +-- +-- betas = +-- { +-- +-- } +-- } +-- } +-- } +--} +function PotentialCSpatialLagRegression(component) + -- Handles with the execution method of a PotentialCSpatialLagRegression component. + -- @arg event A representation of a time instant when the simulation engine must run. + -- @arg luccMEModel A LuccME model. + -- @usage --DONTRUN + -- component.run(event, model) + component.run = function(self, event, luccMEModel) + local luTypes = luccMEModel.landUseTypes + local demand = luccMEModel.demand + local regionsNumber = #self.potentialData + + -- Create an internal constant that can be modified during allocation + for rNumber = 1, regionsNumber, 1 do + for i, luData in pairs(self.potentialData[rNumber]) do + if (luData.const == nil) then + luData.const = 0 + end + + if (luData.minReg == nil) then + luData.minReg = 0 + end + + if (luData.maxReg == nil) then + luData.maxReg = 1 + end + + luData.newconst = luData.const + luData.newminReg = luData.minReg + luData.newmaxReg = luData.maxReg + end + + if (self.constChange == nil) then + self.constChange = 0.1 -- original clue value + end + + if (event:getTime() > luccMEModel.startTime) then + self:adaptRegressionConstants(demand, rNumber) + end + + for i = 1, #luTypes, 1 do + self:computePotential(luccMEModel, rNumber, i) + end + end + end -- function run + + -- Handles with the verify method of a PotentialCSpatialLagRegression component. + -- @arg event A representation of a time instant when the simulation engine must run. + -- @arg luccMEModel A LuccME model. + -- @usage --DONTRUN + -- component.verify(event, self) + component.verify = function(self, event, luccMEModel) + print("Verifying Potential parameters") + local cs = luccMEModel.cs + + if (self.regionAttr == nil) then + self.regionAttr = "region" + end + + forEachCell(cs, function(cell) + cell["alternate_model"] = 0 + + if (cell[self.regionAttr] == nil) then + cell["region"] = 1 + else + cell["region"] = cell[self.regionAttr] + end + end + ) + + -- check potentialData + if (self.potentialData == nil) then + error("potentialData is missing", 2) + end + + local regionsNumber = #self.potentialData + + -- check number of Regions + if (regionsNumber == nil or regionsNumber == 0) then + error("The model must have at least One region") + else + for i = 1, regionsNumber, 1 do + local regressionNumber = #self.potentialData[i] + local lutNumber = #luccMEModel.landUseTypes + + -- check the number of regressions + if (regressionNumber ~= lutNumber) then + error("Invalid number of regressions on Region number "..i.." . Regressions: "..regressionNumber.." LandUseTypes: "..lutNumber) + end + + for j = 1, regressionNumber, 1 do + -- check isLog variable + if(self.potentialData[i][j].isLog == nil) then + error("isLog variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + + -- check minReg variable + if(self.potentialData[i][j].minReg == nil) then + error("minReg variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + + -- check maxReg variable + if(self.potentialData[i][j].maxReg == nil) then + error("maxReg variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + + -- check ro variable + if(self.potentialData[i][j].ro == nil) then + error("ro variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + + -- check const variable + if (self.potentialData[i][j].const == nil) then + error("const variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + + -- check betas variable + if (self.potentialData[i][j].betas == nil) then + error("minReg variable is missing on Region "..i.." LandUseType: "..luccMEModel.landUseTypes[j], 2) + end + + -- check betas within database + for k, lu in pairs (self.potentialData[i][j].betas) do + if (luccMEModel.cs.cells[1][k] == nil) then + error("Beta "..k.." on Region "..i.." LandUseType "..luccMEModel.landUseTypes[j].." not found within database", 2) + end + end + end -- for j + end -- for i + end -- else + + local filename = self.filename + + if (filename ~= nil) then + loadGALNeighborhood(filename) + else + cs:createNeighborhood() + end + end -- verify + + -- Handles with the modify method of a PotentialCSpatialLagRegression component. + -- @arg luccMEModel A LuccME model. + -- @arg rNumber The potential region number. + -- @arg luIndex A land use index (an specific luIndex of a list of possible land uses). + -- @arg direction The direction for the regression. + -- @usage --DONTRUN + -- component.modify(luccMEModel, j, i, luDirect) + component.modify = function(self, luccMEModel, rNumber, luIndex, direction) + local cs = luccMEModel.cs + local luData = self.potentialData[rNumber][luIndex] + + if luData.newconst == nil then + luData.newconst = 0 + end + + if (luData.isLog) then + local const_unlog = (10 ^ luData.newconst) + self.constChange * direction + + if (const_unlog ~= 0) then + luData.newconst = math.log(const_unlog, 10) + end + else + luData.newconst = luData.newconst + self.constChange * direction + end + + self:computePotential (luccMEModel, rNumber, luIndex) + end -- function modifyPotential + + -- Handles with the constants regression method of a PotentialCSpatialLagRegression component. + -- @arg demand A demand to calculate the potential. + -- @arg rNumber The potential region number. + -- @usage --DONTRUN + -- component.adaptRegressionConstants(demand, rNumbers) + component.adaptRegressionConstants = function(self, demand, rNumber) + for i, luData in pairs (self.potentialData[rNumber]) do + local currentDemand = demand:getCurrentLuDemand(i) + local previousDemand = demand:getPreviousLuDemand(i) + local plus = 0.01 * ((currentDemand - previousDemand) / previousDemand) + + luData.newconst = luData.const + + if (luData.isLog) then + local const_unlog = (10 ^ luData.newconst) + plus + + if (const_unlog ~= 0) then + luData.newconst = math.log (const_unlog, 10) + end + else + luData.newconst = luData.newconst + plus + end + + luData.newminReg = luData.newminReg + plus + luData.newmaxReg = luData.newmaxReg + plus + luData.const = luData.newconst + end + end -- function adaptRegressionConstants + + -- Modify potencial for an protected area. + -- @arg complementarLU Land use name. + -- @arg attrProtection The protetion attribute name. + -- @arg rate A rate for potencial multiplier. + -- @arg event A representation of a time instant when the simulation engine must run. + -- @arg luccMEModel A LuccME model. + -- @usage --DONTRUN + -- component.modifyDriver(self.complementarLU, self.attrProtection, 0.5, event, luccMEModel) + component.modifyDriver = function(self, complementarLU, attrProtection, rate, event, luccMEModel) + local regionsNumber = #luccMEModel.potential.potentialData + local luTypes = luccMEModel.landUseTypes + local luIndex = 1 + + for i, complementarLU in pairs (luTypes) do + if (complementarLU == luTypes[i]) then + luIndex = i + break + end + end + + for i = 1, regionsNumber, 1 do + local regressionNumber = #luccMEModel.potential.potentialData[i] + + for j = 1, regressionNumber, 1 do + if (luccMEModel.potential.potentialData[i][j].betas[attrProtection] ~= nil) then + luccMEModel.potential.potentialData[i][j].betas[attrProtection] = luccMEModel.potential.potentialData[i][j].betas[attrProtection] * rate + end + end + + self:computePotential(luccMEModel, i, luIndex) + end + end + + -- Handles with the compute potential method of a PotentialCSpatialLagRegression component. + -- @arg luccMEModel A LuccME model. + -- @arg rNumber The pontencial region number. + -- @arg luIndex A land use index (an specific luIndex of a list of possible land uses). + -- @usage --DONTRUN + -- component.computePotential(luccMEModel, rNumber, luIndex) + component.computePotential = function(self, luccMEModel, rNumber, luIndex) + local cs = luccMEModel.cs + local luTypes = luccMEModel.landUseTypes + local lu = luTypes[luIndex] + local luData = self.potentialData[rNumber][luIndex] + local pot = lu.."_pot" + local reg = lu.."_reg" + local activeRegionNumber = 0 + + for k,cell in pairs (cs.cells) do + if (cell.region == rNumber) then + activeRegionNumber = rNumber + + local regressionX = 0 + + for var, beta in pairs (luData.betas) do + regressionX = regressionX + beta * cell[var] + end + + --regresDrivers = regressionX + local regresY = 0 + local neighSum = 0 + local count = 0 + local neighY = 0 + local Y = 0 + + forEachNeighbor(cell, function(neigh, _, cell) + Y = cell.past[lu] + neighY = neigh.past[lu] + + if (cell[luccMEModel.landUseNoData] ~= 1) then + Y = Y / (1 - cell[luccMEModel.landUseNoData]) + end + + if (neigh[luccMEModel.landUseNoData] ~= 1) then + neighY = neighY / (1 - neigh[luccMEModel.landUseNoData]) + end + + if (neigh[luccMEModel.landUseNoData] < 1) then + count = count + 1 + neighSum = neighSum + neighY + end + end + ) + + if (count > 0) then + regresY = (Y + neighSum) / (count + 1) + else + regresY = Y + end + + local oldRegress = regresY + + if (luData.isLog) then -- if the land use is log transformed + regresY = math.log(regresY + 0.0001, 10) -- ANAP + end + + + regresY = regresY*luData.ro + + local regression = luData.newconst + regressionX + regresY + local regressionLimit = luData.const+ regressionX + regresY + + if (luData.isLog) then -- if the land use is log transformed + regression = (10 ^ regression) - 0.0001 + regressionLimit = (10 ^ regressionLimit) - 0.0001 + end + + local oldReg = regressionLimit + + if (regressionLimit > luData.maxReg) then + regression = 1 + end + + if (regressionLimit < luData.minReg) then + regression = 0 + end + + --if (regression > 1) then -- ANAP + -- regression = 1 + --end + + --if (regression < 0) then + -- regression = 0 + --end + + regression = regression * (1 - cell[luccMEModel.landUseNoData]) + cell[reg] = regression + cell[pot] = regression - cell.past[lu] + end -- if region + end -- for k + + if (activeRegionNumber == 0) then + error("Region ".. rNumber.." is not set into database.") + end + end -- function computePotential + + collectgarbage("collect") + return component +end -- PotentialCSpatialLagRegression diff --git a/tests/test_executor_saturation.py b/tests/test_executor_saturation.py new file mode 100644 index 0000000..b1609fc --- /dev/null +++ b/tests/test_executor_saturation.py @@ -0,0 +1,83 @@ +""" +LuccSaturationExecutor through the local CLI, with the example TOML +(examples/dissmodel-configs/lucc_saturation.toml = lab03), against the lab03 +golden: the path a model like LuccME-BR is run by. + +csAC goes in as a GeoTIFF with named bands (land uses, drivers, mask, the +cells' order), written with dissmodel's save_geotiff; the CLI writes the last +step and `save_steps`; both are compared with the golden. +""" +from __future__ import annotations + +import subprocess +import sys +from pathlib import Path + +import numpy as np +import pytest +from _lab03_helpers import LAND_USES, Lab +from dissmodel.geo.raster.backend import RasterBackend +from dissmodel.io import load_dataset +from dissmodel.io.raster import save_geotiff +from rasterio.transform import from_origin + +ROOT = Path(__file__).resolve().parent.parent +TOML = ROOT / "examples" / "dissmodel-configs" / "lucc_saturation.toml" +DEMAND_CSV = ROOT / "data" / "input" / "examples_demand_lab1.csv" # lab03's demand table, same values +TOL = 1e-9 + + +@pytest.fixture(scope="module") +def run_cli(tmp_path_factory): + lab = Lab("lab03") + tmp = tmp_path_factory.mktemp("saturation") + backend = RasterBackend(shape=lab.shape) + for lu in LAND_USES: + backend.set(lu, lab.initial(lu)) + for name, arr in lab.drivers(2008).items(): + backend.set(name, arr) + backend.set("mask", lab.valid.astype(np.float64)) + backend.set("order", lab.grid(np.arange(len(lab.cells)))) + names = sorted(backend.arrays) + transform = from_origin(0, lab.shape[0] * 5000.0, 5000.0, 5000.0) # any georeference: 5 km cells + cellspace = tmp / "csAC.tif" + save_geotiff( + (backend, {"crs": "EPSG:5880", "transform": transform}), + str(cellspace), + band_spec=[(n, "float64", -9999.0) for n in names], + ) + subprocess.run( + [sys.executable, "-m", "disslucc.executors.saturation", "run", "--toml", str(TOML), + "--input", str(cellspace), "--param", f"demand_csv={DEMAND_CSV}", "--output", str(tmp / "lab03.tif")], + check=True, capture_output=True, text=True, cwd=ROOT, + ) + # the CLI adds the experiment id to the name: lab03_.tif, lab03__step3.tif + outputs = {p.name: p for p in tmp.glob("lab03_*.tif")} + final = next(p for n, p in outputs.items() if "_step" not in n) + return lab, {6: final, 3: final.with_name(f"{final.stem}_step3.tif")} + + +def read(path: Path) -> dict[str, np.ndarray]: + (backend, _), _ = load_dataset(str(path), fmt="raster") + return backend.arrays + + +@pytest.mark.parametrize("step", [3, 6]) +def test_cli_run_matches_golden(run_cli, step): + lab, paths = run_cli + arrays = read(paths[step]) + ref = lab.golden_year(2008 + step) + r = ref.index.get_level_values("row").values + c = ref.index.get_level_values("col").values + for lu in LAND_USES: + err = np.abs(arrays[lu][r, c] - ref[f"{lu}_out"].values).max() + assert err < TOL, f"{2008 + step} {lu}: {err:.3e}" + + +def test_demand_csv_is_lab03s(): + """The example reuses examples_demand_lab1.csv: it holds lab03's demand table.""" + from _lab03_helpers import DEMAND + + from disslucc.components.demand import load_demand_csv + + assert load_demand_csv(DEMAND_CSV.read_text(), LAND_USES) == DEMAND diff --git a/tests/test_lua_differential.py b/tests/test_lua_differential.py new file mode 100644 index 0000000..ad7603a --- /dev/null +++ b/tests/test_lua_differential.py @@ -0,0 +1,312 @@ +""" +The ports against the original LuccME Lua code, run as it is (lupa), on +synthetic cases built to reach what the lab03/lab06 goldens never do: + +- ``correctCellChange`` — in those labs no cell ever needs correcting; +- the saturation branch of ``computeChange`` and ``updateAllocationParameters`` + — those labs set ``changeLimiarValue = 1``; +- ``computePotential`` with log-transformed classes, isolated cells and cells + whose neighbours are all no-data — csAC has none. + +tests/lua/ holds the two Lua files unchanged (LuccME 6244dd4, luccme/lua/ of +terrame-docker v0.1.1) — provisional, see docs/decisions.md. The TerraME +functions they call are stubbed below: ``forEachCell``, ``forEachNeighbor`` (the +neighbourhoods are built here, as TerraME's ``createNeighborhood`` would: Moore +without the cell for the potential, 3 × 3 with the cell for "11x11"), ``belong``. +""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +import pytest + +from disslucc.components.allocation.saturation import ( + compute_change, + correct_cell_change, + saturation_indicator, +) +from disslucc.components.potential.spatial_lag import spatial_lag_regression +from disslucc.schemas import SaturationAllocationSpec, SpatialLagRegressionSpec + +lupa = pytest.importorskip("lupa") + +LUA = Path(__file__).resolve().parent / "lua" +TOL = 1e-12 +STUBS = """ +print = function() end +function forEachCell(cs, f) for _, c in ipairs(cs.cells) do f(c) end end +function forEachNeighbor(cell, name, f) + if type(name) == "function" then f = name; name = "1" end + for _, n in ipairs(cell.neighborhoods[name]) do f(n, 1, cell) end +end +function belong(v, t) for _, x in ipairs(t) do if x == v then return true end end return false end +function makeDemand(dirs) + local d = {dirs = dirs} + function d:getCurrentLuDirection(i) return self.dirs[i] end + return d +end +function makeEvent(t) local e = {t = t}; function e:getTime() return self.t end; return e end +""" + + +@pytest.fixture(scope="module") +def lua(): + rt = lupa.LuaRuntime(unpack_returned_tuples=True) + rt.execute(STUBS) + for name in ("AllocationCClueLikeSaturation.lua", "PotentialCSpatialLagRegression.lua"): + rt.execute((LUA / name).read_text()) + return rt + + +def table(lua, obj): + return lua.table_from(obj, recursive=True) + + +def lua_allocation(lua, specs: list[SaturationAllocationSpec], **extra): + data = [ + { + "static": s.static, + "minValue": s.min_value, + "maxValue": s.max_value, + "minChange": s.min_change, + "maxChange": s.max_change, + "changeLimiarValue": s.change_limiar_value, + "maxChangeAboveLimiar": s.max_change_above_limiar, + } + for s in specs + ] + return lua.globals().AllocationCClueLikeSaturation(table(lua, {"allocationData": [data], **extra})) + + +class Grid: + """A small raster with holes; the cells in row-major order are LuccME's cells.""" + + def __init__(self, rng, shape=(24, 30), holes=0.15, isolated=4): + self.shape = shape + self.valid = rng.random(shape) > holes + # a few isolated cells: a valid cell with every neighbour removed + for _ in range(isolated): + r, c = rng.integers(2, shape[0] - 2), rng.integers(2, shape[1] - 2) + self.valid[r - 1 : r + 2, c - 1 : c + 2] = False + self.valid[r, c] = True + self.rows, self.cols = np.nonzero(self.valid) + + def neighbourhoods(self, lua, cells): + index = {(r, c): k for k, (r, c) in enumerate(zip(self.rows, self.cols, strict=True))} + for k, (r, c) in enumerate(zip(self.rows, self.cols, strict=True)): + moore, window = [], [] + for dr in (-1, 0, 1): + for dc in (-1, 0, 1): + j = index.get((r + dr, c + dc)) + if j is None: + continue + window.append(cells[j]) + if (dr, dc) != (0, 0): + moore.append(cells[j]) + cells[k].neighborhoods = table(lua, {"1": moore, "11x11": window}) + + +def random_specs(rng, n) -> list[SaturationAllocationSpec]: + specs = [] + for i in range(n): + lo = float(rng.choice([0.0, 0.0, 0.05, 0.2])) + hi = float(rng.choice([1.0, 1.0, 0.8, 0.6])) + specs.append( + SaturationAllocationSpec( + static=int(rng.choice([-1, -1, 0, 1])) if i else -1, + min_value=lo, + max_value=hi, + min_change=float(rng.choice([0.0, 0.01])), + max_change=float(rng.choice([1.0, 0.3, 0.1])), + change_limiar_value=float(rng.uniform(0.2, 0.8)), + max_change_above_limiar=float(rng.uniform(0.0, 0.1)), + ) + ) + return specs + + +# ── correctCellChange ────────────────────────────────────────────────────── + + +@pytest.mark.parametrize("seed", range(12)) +def test_correct_cell_change_matches_lua(lua, seed): + rng = np.random.default_rng(seed) + n_cells, n_lu = 300, int(rng.integers(3, 5)) + lus = [f"c{i}" for i in range(n_lu)] + specs = random_specs(rng, n_lu) + directions = [int(d) for d in rng.choice([-1, 0, 1], size=n_lu)] + + past = rng.dirichlet(np.ones(n_lu), size=n_cells) + step = rng.normal(0, 0.08, size=(n_cells, n_lu)) + same_way = rng.random(n_cells) < 0.3 # every class moving the same way: the "shift" branch + step[same_way] = np.abs(step[same_way]) * rng.choice([-1, 1], size=(same_way.sum(), 1)) + values = np.clip(past + step, 0, 1) + + cells = [] + for k in range(n_cells): + cell = {lu: float(values[k, i]) for i, lu in enumerate(lus)} + cell["past"] = {lu: float(past[k, i]) for i, lu in enumerate(lus)} + cell["regionAloc"] = 1 + cells.append(table(lua, cell)) + model = table(lua, {"cs": {"cells": cells}, "landUseTypes": lus}) + model.demand = lua.globals().makeDemand(table(lua, directions)) + comp = lua_allocation(lua, specs) + comp.correctCellChange(comp, model, 1) + theirs = np.array([[cells[k][lu] for lu in lus] for k in range(n_cells)]) + + stats: dict = {} + ours = correct_cell_change(values, past, specs, directions, stats) + assert np.abs(ours - theirs).max() < TOL, stats + assert stats["need"] > 0 + + +def test_correct_cell_change_cases_reach_every_branch(lua): + """The random cases above are only worth something if they reach the branches.""" + total: dict = {} + for seed in range(12): + rng = np.random.default_rng(seed) + n_lu = int(rng.integers(3, 5)) + specs = random_specs(rng, n_lu) + past = rng.dirichlet(np.ones(n_lu), size=300) + step = rng.normal(0, 0.08, size=(300, n_lu)) + same_way = rng.random(300) < 0.3 + step[same_way] = np.abs(step[same_way]) * rng.choice([-1, 1], size=(same_way.sum(), 1)) + values = np.clip(past + step, 0, 1) + correct_cell_change(values, past, specs, [int(d) for d in rng.choice([-1, 0, 1], size=n_lu)], total) + assert total["need"] > 1000 + assert total["shift"] > 100 + assert total["backp_lowered"] > 100 + assert total["several_rounds"] > 10 + + +# ── computeChange with saturation ────────────────────────────────────────── + + +@pytest.mark.parametrize("seed", range(8)) +def test_compute_change_matches_lua(lua, seed): + rng = np.random.default_rng(100 + seed) + n_cells, n_lu = 400, 3 + lus = ["f", "d", "o"] + specs = random_specs(rng, n_lu) + directions = [int(d) for d in rng.choice([-1, 1], size=n_lu)] + elasticity = [float(e) for e in rng.uniform(0.05, 1.5, size=n_lu)] + past = rng.dirichlet(np.ones(n_lu), size=n_cells) + pot = rng.normal(0, 0.3, size=(n_cells, n_lu)) + pot[rng.random((n_cells, n_lu)) < 0.05] = 0.0 + saturation = rng.random(n_cells) + + cells = [] + for k in range(n_cells): + cell = {f"{lu}_pot": float(pot[k, i]) for i, lu in enumerate(lus)} + cell.update({lu: float(past[k, i]) for i, lu in enumerate(lus)}) + cell["past"] = {lu: float(past[k, i]) for i, lu in enumerate(lus)} + cell["regionAloc"] = 1 + cell["sat"] = float(saturation[k]) + cells.append(table(lua, cell)) + model = table(lua, {"cs": {"cells": cells}, "landUseTypes": lus}) + model.demand = lua.globals().makeDemand(table(lua, directions)) + comp = lua_allocation(lua, specs, saturationIndicator="sat") + comp.elasticity = table(lua, elasticity) + comp.computeChange(comp, model, 1) + + saturated = 0 + for i, lu in enumerate(lus): + theirs = np.array([cells[k][lu] for k in range(n_cells)]) + ours = compute_change(past[:, i], pot[:, i], elasticity[i], specs[i], directions[i], saturation) + assert np.abs(ours - theirs).max() < TOL, (lu, specs[i]) + saturated += int((saturation > specs[i].change_limiar_value).sum()) + assert saturated > 0 + + +# ── updateAllocationParameters (saturation indicator) ────────────────────── + + +@pytest.mark.parametrize("seed", range(6)) +def test_saturation_indicator_matches_lua(lua, seed): + rng = np.random.default_rng(200 + seed) + grid = Grid(rng) + shape = grid.shape + comp_share = rng.random(shape) + no_data = np.where(rng.random(shape) < 0.1, 1.0, rng.random(shape) * 0.3) + protection = np.where(rng.random(shape) < 0.2, rng.random(shape) * 0.8, 0.0) + + cells = [] + for r, c in zip(grid.rows, grid.cols, strict=True): + cells.append( + table( + lua, + {"f": float(comp_share[r, c]), "nd": float(no_data[r, c]), "prot": float(protection[r, c])}, + ) + ) + grid.neighbourhoods(lua, cells) + model = table(lua, {"cs": {"cells": cells}, "landUseNoData": "nd"}) + comp = lua_allocation( + lua, + [SaturationAllocationSpec()], + saturationIndicator="sat", + attrProtection="prot", + complementarLU="f", + ) + comp.updateAllocationParameters(comp, lua.globals().makeEvent(2000), model) + + theirs = np.array([cell["sat"] for cell in cells]) + ours = saturation_indicator(comp_share, grid.valid, no_data, protection)[grid.rows, grid.cols] + assert np.abs(ours - theirs).max() < TOL + + +# ── computePotential ─────────────────────────────────────────────────────── + + +@pytest.mark.parametrize("is_log", [False, True]) +@pytest.mark.parametrize("seed", range(4)) +def test_spatial_lag_potential_matches_lua(lua, seed, is_log): + rng = np.random.default_rng(300 + seed) + grid = Grid(rng) + shape = grid.shape + share = rng.random(shape) + no_data = np.where(rng.random(shape) < 0.25, 1.0, rng.random(shape) * 0.4) # many all-no-data cells + drivers = {"x1": rng.random(shape), "x2": rng.normal(0, 1, shape)} + spec = SpatialLagRegressionSpec( + const=float(rng.normal(0, 0.2)), + ro=float(rng.uniform(0.5, 1.0)), + betas={"x1": float(rng.normal(0, 0.5)), "x2": float(rng.normal(0, 0.1))}, + is_log=is_log, + min_reg=float(rng.choice([0.0, 0.1])), + max_reg=float(rng.choice([1.0, 0.9])), + ) + newconst = spec.const + 0.1 * int(rng.integers(-2, 3)) # as after the allocation's modify + + cells = [] + for r, c in zip(grid.rows, grid.cols, strict=True): + cell = {"region": 1, "nd": float(no_data[r, c]), "past": {"u": float(share[r, c])}} + cell.update({k: float(v[r, c]) for k, v in drivers.items()}) + cells.append(table(lua, cell)) + grid.neighbourhoods(lua, cells) + model = table(lua, {"cs": {"cells": cells}, "landUseTypes": ["u"], "landUseNoData": "nd"}) + data = { + "isLog": is_log, + "const": spec.const, + "minReg": spec.min_reg, + "maxReg": spec.max_reg, + "ro": spec.ro, + "betas": spec.betas, + } + comp = lua.globals().PotentialCSpatialLagRegression(table(lua, {"potentialData": [[data]]})) + comp.potentialData[1][1].newconst = newconst + comp.computePotential(comp, model, 1, 1) + + theirs = np.array([cell["u_pot"] for cell in cells]) + _, pot = spatial_lag_regression(spec, share, drivers, grid.valid, no_data, newconst) + ours = pot[grid.rows, grid.cols] + assert np.abs(ours - theirs).max() < 1e-10 + + # the cases do contain what csAC lacks + moore_count = sum( + np.roll(np.roll(grid.valid, dr, 0), dc, 1) + for dr in (-1, 0, 1) + for dc in (-1, 0, 1) + if (dr, dc) != (0, 0) + ) + assert ((moore_count == 0) & grid.valid).any() diff --git a/tests/test_saturation_golden.py b/tests/test_saturation_golden.py new file mode 100644 index 0000000..881218b --- /dev/null +++ b/tests/test_saturation_golden.py @@ -0,0 +1,62 @@ +""" +AllocationClueLikeSaturation + PotentialSpatialLagRegression against the LuccME +goldens of lab03 and lab06: the whole model, run from 2008 on its own (no replay). + +Checked every year, 2008–2014: +- every cell's land use (_out) and potential (_pot) against the golden; +- the maximum error against the demand, and the iteration count, against the + TerraME log (terrame.log: "Maximum error: …", "Number of iterations: …"). + +The saturation branch is not reached in these labs (change_limiar_value = 1), +so it is not constrained here; everything else in the allocation is. +""" + +from __future__ import annotations + +import pytest +from _lab03_helpers import ALLOCATION, YEARS, Lab, run, terrame_log, worst_per_column + +TOL = 1e-9 + + +@pytest.fixture(scope="module", params=["lab03", "lab06"]) +def result(request): + lab = Lab(request.param) + snaps, allocation = run(lab) + return lab, snaps, allocation + + +def test_land_use_and_potential_match_golden_every_year(result): + lab, snaps, _ = result + worst = worst_per_column(lab, snaps) + bad = {col: w for col, w in worst.items() if w[1] >= TOL} + assert not bad, f"{lab.name}: {bad}" + + +def test_iterations_and_error_match_terrame_log(result): + lab, _, allocation = result + log = terrame_log(lab.name) + assert sorted(log) == YEARS + assert allocation.iterations_per_step == [log[y][0] for y in YEARS] + for year, ours in zip(YEARS, allocation.max_error_per_step, strict=True): + assert ours == pytest.approx(log[year][1], rel=1e-9), year + + +def test_what_the_goldens_do_not_reach(monkeypatch): + """In lab03 no cell ever needs correctCellChange (the classes always sum to + 1 ± 0.005 before it) and no cell is saturated (change_limiar_value = 1): the + goldens say nothing about those parts, and the cells' visiting order has no + effect. They are checked against the Lua itself in test_lua_differential.py.""" + import disslucc.components.allocation.saturation as sat + + stats: dict = {} + original = sat.correct_cell_change + monkeypatch.setattr(sat, "correct_cell_change", lambda v, p, s, d: original(v, p, s, d, stats)) + lab = Lab("lab03") + run(lab) + assert stats["need"] == 0 + assert all(s.change_limiar_value >= 1 for s in ALLOCATION) + + monkeypatch.setattr(sat, "correct_cell_change", original) + snaps, _ = run(lab, order="shuffled") + assert max(err for _, err in worst_per_column(lab, snaps).values()) < TOL diff --git a/tests/test_spatial_lag_golden.py b/tests/test_spatial_lag_golden.py new file mode 100644 index 0000000..26d504e --- /dev/null +++ b/tests/test_spatial_lag_golden.py @@ -0,0 +1,125 @@ +""" +PotentialSpatialLagRegression against the LuccME goldens of lab03 and lab06. + +Goldens (benchmark/goldens/, from LambdaGeo/terrame-docker v0.1.1): the state of every +cell at the end of every year, 2008–2014, with _out and _pot. Both labs +use DemandPreComputedValues + PotentialCSpatialLagRegression + +AllocationCClueLikeSaturation; lab06 also replaces `ti` from csAC_2009 in 2009. + +To test the potential alone, each year starts from the golden's own land use of +the year before (what LuccME's `cell.past` holds) — so an error cannot come from +the allocation. The allocation can still shift `newconst` by ±0.1 during a year +(`modify`); in these two labs it never does, so the potential must match exactly +(the golden keeps 12 decimals: tolerance 1e-10). +""" + +from __future__ import annotations + +import numpy as np +import pytest +from _lab03_helpers import DEMAND, LAND_USES, NO_DATA, YEARS, Lab, max_errors, specs +from dissmodel.core import Environment, Model +from dissmodel.geo.raster.backend import RasterBackend + +from disslucc import DemandPreComputedValues, PotentialSpatialLagRegression, SpatialLagRegressionSpec +from disslucc.components.potential.spatial_lag import spatial_lag_regression + +TOL = 1e-10 + + +@pytest.fixture(scope="module", params=["lab03", "lab06"]) +def lab(request) -> Lab: + return Lab(request.param) + + +def test_golden_is_the_expected_run(lab): + assert lab.manifest["status"] == "ok" + assert lab.manifest["years"] == [YEARS[0], YEARS[-1]] + assert lab.manifest["n_cells"] == len(lab.cells) == len(lab.golden_year(2008)) + + +def test_potential_matches_golden_every_year_and_class(lab): + errors = max_errors(lab) + worst = max(errors, key=errors.get) + assert errors[worst] < TOL, f"{lab.name}: worst {worst} = {errors[worst]:.3e}" + + +def test_clipping_is_exercised(lab): + """The golden constrains the clip branch only if some cells hit it: count the + cells whose value changes when the bounds are removed.""" + past, drivers = lab.past(2008), lab.drivers(2008) + clipped = 0 + for lu, spec in zip(LAND_USES, specs(), strict=True): + free = SpatialLagRegressionSpec( + const=spec.const, ro=spec.ro, betas=spec.betas, min_reg=-np.inf, max_reg=np.inf + ) + reg, _ = spatial_lag_regression(spec, past[lu], drivers, lab.valid, past[NO_DATA]) + reg_free, _ = spatial_lag_regression(free, past[lu], drivers, lab.valid, past[NO_DATA]) + clipped += int((np.abs(reg - reg_free)[lab.valid] > 1e-12).sum()) + assert clipped > 0 + + +def test_non_cumulative_adaptation_is_rejected(lab): + errors = max_errors(lab, cumulative=False) + assert max(v for (y, _), v in errors.items() if y >= 2010) > 1e-4 + + +def test_updates_pair_cells_by_position_not_by_id(): + lab = Lab("lab06") + by_id = max_errors(lab, pair_by="id") + assert max(v for (y, _), v in by_id.items() if y >= 2009) > 1e-2 + + +def test_component_matches_golden(lab): + """The dissmodel component (what moves to disslucc), stepped in an Environment. + + A replay model runs before it each year and puts the golden's land use of the + year before into the backend, plus the drivers of that year.""" + backend = RasterBackend(shape=lab.shape) + for lu in LAND_USES: + backend.set(lu, lab.initial(lu)) + for name, arr in lab.drivers(YEARS[0]).items(): + backend.set(name, arr) + backend.set("mask", lab.valid.astype(np.float32)) + + class Replay(Model): + def setup(self): + self.year = YEARS[0] + + def execute(self): + self.year = YEARS[0] + int(self.env.now()) + if self.year > YEARS[0]: + for lu, arr in lab.past(self.year).items(): + backend.set(lu, arr) + backend.set(lu + "_past", arr.copy()) + for name, arr in lab.drivers(self.year).items(): + backend.set(name, arr) + + got: dict[int, dict[str, np.ndarray]] = {} + + class Record(Model): + def execute(self): + year = YEARS[0] + int(self.env.now()) + got[year] = {lu: backend.get(lu + "_pot").copy() for lu in LAND_USES} + + env = Environment(end_time=len(YEARS) - 1) + demand = DemandPreComputedValues(annual_demand=DEMAND, land_use_types=LAND_USES) + Replay() + PotentialSpatialLagRegression( + backend=backend, + potential_data=[specs()], + demand=demand, + land_use_types=LAND_USES, + land_use_no_data=NO_DATA, + ) + Record() + env.run() + + assert sorted(got) == YEARS + for year in YEARS: + ref = lab.golden_year(year) + r = ref.index.get_level_values("row").values + c = ref.index.get_level_values("col").values + for lu in LAND_USES: + err = np.abs(got[year][lu][r, c] - ref[f"{lu}_pot"].values).max() + assert err < TOL, f"{lab.name} {year} {lu}: {err:.3e}"