diff --git a/README.md b/README.md index 5fc0d50..85de34e 100644 --- a/README.md +++ b/README.md @@ -203,7 +203,7 @@ Infrastructure Layer (ModelExecutor) ### Benchmark executor -The `LUCCBenchmarkExecutor` is a meta-executor that runs vector and raster substrates in a single pass and compares both against a TerraME/LUCCME reference result. It generates a Markdown report and scatter plots — the primary tool for validating numerical equivalence before publishing results. +The `LUCCBenchmarkExecutor` is a meta-executor that runs vector and raster substrates in a single pass and compares both against a TerraME/LUCCME reference result. It generates a Markdown report and scatter plots — the primary tool for quantifying agreement with the reference before publishing results. See “Interpreting the numbers” below: the output is an agreement measurement, not a claim of cell-by-cell equivalence. ```bash python -m disslucc_continuous.executors.lucc_benchmark_executor run \ @@ -400,7 +400,7 @@ DisSLUCC-Continuous/ 3. **Reproducibility** — Each experiment records model commit, input checksum, and resolved spec via `ExperimentRecord`. 4. **Two substrates** — Same algorithms available for vector (GeoDataFrame) and raster (RasterBackend/NumPy). 5. **Executor pattern** — Science layer never knows about files or cloud; infrastructure layer never calculates spatial equations. -6. **Benchmark-first validation** — `lucc_benchmark` validates numerical equivalence between substrates and against TerraME/LUCCME reference results before any production use. +6. **Benchmark-first validation** — `lucc_benchmark` quantifies agreement against TerraME/LUCCME reference results, and `tests/test_benchmark_discriminance.py` verifies that the benchmark actually constrains the implementation, before any production use. --- @@ -412,7 +412,7 @@ The primary validation strategy is numerical comparison against TerraME/LUCCME r |---|---| | `Vector_vs_TerraME` | Python vector model vs original LUCCME/TerraME result | | `Raster_vs_TerraME` | Python raster model vs original LUCCME/TerraME result | -| `Vector_vs_Raster` | Internal consistency between substrates | +| `Vector_vs_Raster` | Consistency of array reshaping between substrates (see caveat below) | Run the full test suite (benchmark validation + unit tests): @@ -425,6 +425,56 @@ The benchmark test (`tests/test_benchmark_validation.py`) uses: - **Reference**: `benchmark/data/LUCCME_Lab1_2014.zip` - **Assertion**: MAE and RMSE below `tolerance=0.01` for `Vector_vs_TerraME` and `Raster_vs_TerraME` +`tests/test_benchmark_discriminance.py` complements this by checking that the +benchmark actually constrains the implementation: perturbing the regression +coefficients (halved, doubled, zeroed, sign-flipped) must break the tolerance +criterion, and a do-nothing baseline must fail it too. All of these hold — the +Lab1 scenario is genuinely discriminative. + +### Interpreting the numbers + +| Metric | What it means | +|---|---| +| `MAE` / `RMSE` | Averages. Low values coexist with a sizeable tail outside tolerance — never report them alone. | +| `Match %` | Fraction of cells within `tolerance`. Currently **87.37%** for `Vector_vs_TerraME`, i.e. ~830 cells (12.63%) differ by more than 0.01, with a max error of 0.027. | +| `Quantity` / `Allocation` | Pontius & Millones (2011) decomposition. `Quantity` is error in the total allocated; `Allocation` is error in *where* it was put. `Quantity + Allocation = MAE`. | + +Because of the above, the accurate phrasing for this result is **"agreement with +MAE below 0.01, with 87% of cells within ±0.01 and a maximum error of 0.027"** — +not "numerical equivalence", which would imply cell-by-cell parity. + +> ### ⚠️ `Vector_vs_Raster` is not a spatial check +> +> The exact agreement between substrates is real — they are independent runs, with +> separate `Environment` instances and a raster backend built from actual `row`/`col` +> values. But there are **no neighbourhood operations anywhere** in `src/`: both +> substrates are purely element-wise arithmetic over the same numbers, merely +> reshaped. This row confirms that NumPy reshapes correctly; it says nothing about +> spatial behaviour, and should not be read as a third independent validation. + +> ### Convergence tolerance — read before interpreting `n_steps` sweeps +> +> The original LuccME script declares `maxDifference = 1643` in +> `AllocationCClueLike`, in area units, against a 2014 demand of 21607.38 for +> class `d` — a **7.6% convergence band**. The Python defaults match it exactly. +> +> This band is wide enough that the reference itself stops 1001.45 area units +> short of the demand it declares. That is legitimate convergence slack, not a +> defect. It also means Python and TerraME can each halt at different, equally +> valid points inside the band, which is why the fit varies non-monotonically +> with `n_steps` (best at 4, official at 6). Sweeps over `n_steps` should not be +> read as evidence of temporal misalignment. +> +> Consistently with this, the Pontius decomposition at the official +> configuration is ~90% quantity and ~10% allocation: the model places +> deforestation in nearly the right cells and misses on the total. Both facts +> are asserted in `tests/test_benchmark_discriminance.py`. + +### Coverage limits + +The comparison covers **only the `d` band at the final step**. The `f` and `outros` +classes are never compared against TerraME, and neither are intermediate steps. + --- ## 📚 References diff --git a/src/disslucc_continuous/executors/lucc_benchmark_executor.py b/src/disslucc_continuous/executors/lucc_benchmark_executor.py index bc81e2f..1f88867 100644 --- a/src/disslucc_continuous/executors/lucc_benchmark_executor.py +++ b/src/disslucc_continuous/executors/lucc_benchmark_executor.py @@ -240,12 +240,35 @@ def save(self, result: dict, record: ExperimentRecord) -> ExperimentRecord: # ── helpers ─────────────────────────────────────────────────────────────────── def _metrics(a: np.ndarray, b: np.ndarray, tol: float) -> dict: + """Agreement metrics between two continuous fields. + + Includes the Pontius & Millones (2011) decomposition adapted to the + continuous case, splitting the error into two components: + + * ``quantity_disagreement`` — difference in totals. Answers "did the model + allocate the right *amount*?". + * ``allocation_disagreement`` — residual error after matching the totals. + Answers "did the model put it in the right *place*?". + + The identity ``quantity + allocation == mae`` holds. A low MAE can mask + either a systematic quantity bias or a spatially wrong allocation, which + is why the decomposition matters. + """ diff = np.abs(a - b) + mae = float(diff.mean()) + + # Pontius & Millones (2011), continuous form. + # |mean(a) - mean(b)| <= mean|a - b|, so allocation is always >= 0. + quantity = float(abs(a.mean() - b.mean())) + allocation = float(mae - quantity) + return { "match_pct": float((diff <= tol).mean() * 100), - "mae": float(diff.mean()), + "mae": mae, "rmse": float(np.sqrt((diff**2).mean())), "max_err": float(diff.max()), + "quantity_disagreement": quantity, + "allocation_disagreement": allocation, "n_cells": len(a), } @@ -322,14 +345,27 @@ def _build_markdown(n_steps: int, tol: float, vec_ms: float, ras_ms: float, metr f"| Vector | {vec_ms:.1f} | 1× |\n", f"| Raster | {ras_ms:.1f} | {vec_ms/ras_ms:.1f}× |\n\n", "## Accuracy — `d`\n\n", - "| Comparison | Match % | MAE | RMSE | Max err | N cells |\n", - "|---|---|---|---|---|---|\n", + "| Comparison | Match % | MAE | Quantity | Allocation | RMSE | Max err | N cells |\n", + "|---|---|---|---|---|---|---|---|\n", ] for label, m in metrics.items(): lines.append( f"| {label.replace('_', ' ')} | {m['match_pct']:.2f}% | {m['mae']:.6f} | " + f"{m['quantity_disagreement']:.6f} | {m['allocation_disagreement']:.6f} | " f"{m['rmse']:.6f} | {m['max_err']:.6f} | {m['n_cells']} |\n" ) + lines.append( + "\n> **Reading these numbers.** `Match %` and `MAE` answer different questions —\n" + "> a low MAE coexists with a sizeable tail outside tolerance, so report both.\n" + "> `Quantity` and `Allocation` are the Pontius & Millones (2011) decomposition:\n" + "> error in *how much* was allocated vs. error in *where*.\n" + "> `Quantity + Allocation = MAE`.\n" + ">\n" + "> **`Vector vs Raster` is not a spatial check.** Neither substrate uses\n" + "> neighbourhood operations — both are element-wise arithmetic over the same\n" + "> values, merely reshaped. Exact agreement here confirms consistent reshaping,\n" + "> not spatial correctness.\n" + ) return "".join(lines) diff --git a/tests/test_benchmark_discriminance.py b/tests/test_benchmark_discriminance.py new file mode 100644 index 0000000..87bbbf3 --- /dev/null +++ b/tests/test_benchmark_discriminance.py @@ -0,0 +1,232 @@ +""" +tests/test_benchmark_discriminance.py +===================================== +Measures the **discriminative power** of the Lab1 benchmark. + +The question here is not "does the model reproduce TerraME?" but "can the +benchmark tell a correct implementation from an incorrect one?". Reproducing the +reference is necessary but not sufficient: if the test still passes with wrong +coefficients, it does not constrain the implementation. + +Known state (2026-07-27) +------------------------ +**Lab1 is discriminative.** Every perturbation tested on the regression +coefficients breaks the tolerance criterion, including zeroing the drivers and +halving the betas. A do-nothing baseline fails too. Unlike the Lab6 scenario in +``disslucc-discrete``, there is no trivial shortcut. + +On the convergence tolerance +---------------------------- +The original LuccME script (``lab1.lua`` / ``lab1_submodel.lua``) declares +``maxDifference = 1643`` in ``AllocationCClueLike`` — in area units, against a +2014 demand of 21607.38 for class ``d``. That is a **7.6% convergence band**, +and the Python defaults match it exactly. + +This band is wide enough that the reference itself stops 1001.45 area units +short of the demand it declares, which is legitimate rather than a defect. It +also means Python and TerraME can each halt at different, equally valid points +inside the band — which is why the fit varies non-monotonically with +``n_steps``. See ``test_reference_gap_within_luccme_tolerance``. +""" +from __future__ import annotations + +import contextlib +import copy +import io +import pathlib + +import numpy as np +import pytest + +ROOT = pathlib.Path(__file__).parent.parent +INPUT_ZIP = ROOT / "examples" / "data" / "input" / "csAC.zip" +DEMAND_CSV = ROOT / "examples" / "data" / "input" / "examples_demand_lab1.csv" +TERRAME_ZIP = ROOT / "benchmark" / "data" / "LUCCME_Lab1_2014.zip" + +TOLERANCE = 0.01 +N_STEPS = 6 +CELL_AREA = 25.0 + +# Values transcribed from the original LuccME script (lab1_submodel.lua). +LUCCME_MAX_DIFFERENCE = 1643.0 # AllocationCClueLike, area units +LUCCME_DEMAND_D_2014 = 21607.38493 # D1.annualDemand, last row, class "d" + +skip_if_no_data = pytest.mark.skipif( + not INPUT_ZIP.exists() or not TERRAME_ZIP.exists(), + reason="Lab1 data files not found", +) + + +def _run_benchmark(n_steps: int = N_STEPS) -> dict: + """Run the benchmark and return the metrics dictionary.""" + from dissmodel.executor import ExperimentRecord + + from disslucc_continuous.executors import LUCCBenchmarkExecutor + + executor = LUCCBenchmarkExecutor() + record = ExperimentRecord( + model_name="lucc_benchmark", + source={"uri": str(INPUT_ZIP)}, + parameters={ + "demand_csv": str(DEMAND_CSV), + "terrame_reference": str(TERRAME_ZIP), + "n_steps": n_steps, + "tolerance": TOLERANCE, + }, + ) + with contextlib.redirect_stdout(io.StringIO()): + result = executor.run(executor.load(record), record) + return result["metrics"] + + +@contextlib.contextmanager +def _scaled_betas(factor: float): + """Temporarily multiply every regression beta by ``factor``.""" + import disslucc_continuous.executors.lucc_benchmark_executor as B + + original = B.POTENTIAL_DATA + perturbed = copy.deepcopy(original) + for spec in perturbed[0]: + if spec.betas: + spec.betas = {c: v * factor for c, v in spec.betas.items()} + B.POTENTIAL_DATA = perturbed + try: + yield + finally: + B.POTENTIAL_DATA = original + + +# ══════════════════════════════════════════════════════════════════════════════ +# 1. Discriminance — perturbing the coefficients must break the test +# ══════════════════════════════════════════════════════════════════════════════ + +@skip_if_no_data +@pytest.mark.parametrize( + "factor,description", + [ + (0.0, "betas zeroed (no drivers)"), + (0.5, "betas halved"), + (2.0, "betas doubled"), + (-1.0, "betas sign-flipped"), + ], +) +def test_perturbed_coefficients_fail_tolerance(factor, description): + """With wrong coefficients the benchmark must fail. + + If it passed, the tolerance would be too loose to constrain the + implementation — and the reported parity would carry no evidential weight. + """ + with _scaled_betas(factor): + m = _run_benchmark()["Vector_vs_TerraME"] + + assert m["mae"] >= TOLERANCE or m["rmse"] >= TOLERANCE, ( + f"With {description}, the benchmark still passed " + f"(MAE={m['mae']:.6f}, RMSE={m['rmse']:.6f}). " + "The tolerance does not constrain the implementation." + ) + + +@skip_if_no_data +def test_do_nothing_baseline_fails(): + """Keeping `d` at its initial value must fail the tolerance criterion.""" + import geopandas as gpd + + gdf = gpd.read_file(str(INPUT_ZIP)) + ref = gpd.read_file(str(TERRAME_ZIP)) + + d_initial = gdf["d"].values.astype(float) + d_out = ref["d_out"].values.astype(float) + mae = float(np.abs(d_initial - d_out).mean()) + + assert mae >= TOLERANCE, ( + f"The do-nothing baseline has MAE={mae:.6f}, below the {TOLERANCE} " + "tolerance. The benchmark would be satisfied without simulating anything." + ) + + +# ══════════════════════════════════════════════════════════════════════════════ +# 2. Characterisation — what the benchmark does NOT cover +# ══════════════════════════════════════════════════════════════════════════════ + +@skip_if_no_data +def test_match_pct_is_reported_alongside_mae(): + """A low MAE coexists with a sizeable tail outside tolerance. + + This test imposes no floor on `match_pct` — it records that the two numbers + tell different stories, and that reporting MAE alone overstates agreement. + It anchors the wording used in the README. + """ + m = _run_benchmark()["Vector_vs_TerraME"] + + assert m["mae"] < TOLERANCE, "regression: the model no longer meets the criterion" + assert m["match_pct"] < 100.0, ( + "If match_pct reached 100%, the README wording can be strengthened — " + "it currently states explicitly that agreement is not cell-by-cell." + ) + + +@skip_if_no_data +def test_pontius_identity_holds(): + """quantity + allocation must reproduce MAE exactly.""" + for label, m in _run_benchmark().items(): + total = m["quantity_disagreement"] + m["allocation_disagreement"] + assert total == pytest.approx(m["mae"], abs=1e-12), ( + f"{label}: Pontius identity violated " + f"({total:.9f} != {m['mae']:.9f})" + ) + assert m["allocation_disagreement"] >= -1e-12, ( + f"{label}: negative allocation disagreement — impossible" + ) + + +# ══════════════════════════════════════════════════════════════════════════════ +# 3. Convergence tolerance — why the fit varies with n_steps +# ══════════════════════════════════════════════════════════════════════════════ + +@skip_if_no_data +def test_reference_gap_within_luccme_tolerance(): + """The reference falls short of its declared demand — legitimately. + + ``AllocationCClueLike`` in the original script declares + ``maxDifference = 1643`` area units. The reference stops ~1001 units short + of the 2014 demand for ``d``, which is inside that band. This is correct + behaviour, not a defect in the reference. + + It also explains why the fit varies non-monotonically with ``n_steps``: + a 7.6% convergence band lets Python and TerraME halt at different, equally + valid points. Any comparison against this reference inherits that slack, + so ``n_steps`` sweeps should not be read as evidence of misalignment. + """ + import geopandas as gpd + + ref = gpd.read_file(str(TERRAME_ZIP)) + ref_area = float(ref["d_out"].values.astype(float).sum()) * CELL_AREA + gap = abs(ref_area - LUCCME_DEMAND_D_2014) + + assert gap < LUCCME_MAX_DIFFERENCE, ( + f"The reference misses its declared demand by {gap:.2f} area units, " + f"beyond the maxDifference={LUCCME_MAX_DIFFERENCE} declared in the " + "LuccME script. That would indicate a genuine problem in the reference " + "rather than ordinary convergence slack." + ) + + +@skip_if_no_data +def test_quantity_error_dominates_allocation_error(): + """Most of the residual is *how much*, not *where*. + + At the official configuration the error decomposes into roughly 90% + quantity and 10% allocation: the model places deforestation in nearly the + right cells and misses on the total. That is the expected signature of a + wide convergence band, and it is the reading that plain MAE hides. + + If this ever inverts — allocation dominating — the cause is spatial and + warrants investigation, because the convergence band does not explain it. + """ + m = _run_benchmark()["Vector_vs_TerraME"] + + assert m["quantity_disagreement"] > m["allocation_disagreement"], ( + f"Allocation disagreement ({m['allocation_disagreement']:.6f}) now " + f"exceeds quantity ({m['quantity_disagreement']:.6f}). The residual is " + "spatial and is no longer explained by the convergence tolerance." + )