From 557aaa5fc56ad79b41a8d63f897f376e13e3503a Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Wed, 16 Sep 2026 17:50:40 -0300 Subject: [PATCH 01/12] Add the calibration objective and the observed series readers rubem.calibration.objective computes the Nash-Sutcliffe efficiency over the pairs valid in both series (NaN, -9999 and the PCRaster missing value masked, a zero-variance guard), the objective 1000 (100 (1 - NSE))^2 and the alignment of a simulated and an observed series on their shared time steps and station ids after a spin-up window; read_series reads the CSV layout of the model's time series tables and PCRaster TSS files. --- rubem/calibration/__init__.py | 12 + rubem/calibration/objective.py | 372 +++++++++++++++++++++ tests/unit/calibration/__init__.py | 0 tests/unit/calibration/test_objective.py | 396 +++++++++++++++++++++++ 4 files changed, 780 insertions(+) create mode 100644 rubem/calibration/__init__.py create mode 100644 rubem/calibration/objective.py create mode 100644 tests/unit/calibration/__init__.py create mode 100644 tests/unit/calibration/test_objective.py diff --git a/rubem/calibration/__init__.py b/rubem/calibration/__init__.py new file mode 100644 index 00000000..35a716c9 --- /dev/null +++ b/rubem/calibration/__init__.py @@ -0,0 +1,12 @@ +"""Automatic calibration of the model parameters against observed streamflow. + +The package holds the objective function (Nash-Sutcliffe efficiency over the +sample stations), the decision vector of the free calibration parameters, the +worker that evaluates one candidate, and the runner that drives SciPy's +differential evolution. + +Importing this package does not import SciPy: the optional dependency is +imported inside the functions that need it, so that reading a series or +computing an efficiency works in an installation without the +``rubem[calibration]`` extra. +""" diff --git a/rubem/calibration/objective.py b/rubem/calibration/objective.py new file mode 100644 index 00000000..7839851e --- /dev/null +++ b/rubem/calibration/objective.py @@ -0,0 +1,372 @@ +"""The calibration objective: Nash-Sutcliffe efficiency over the sample stations. + +The module reads the observed and the simulated time series, masks the values +that stand for a missing observation, computes the Nash-Sutcliffe efficiency +(NSE) of each station and turns the mean efficiency into the objective the +optimizer minimizes. + +Nothing here imports SciPy: the efficiency and the objective are evaluated with +NumPy alone. +""" + +import csv +import logging +from dataclasses import dataclass +from pathlib import Path + +import numpy as np + +from .._paths import PathInput, as_path + +logger = logging.getLogger(__name__) + +MISSING_VALUES = (-9999.0,) +"""Values written by the model and by the usual observation tables for a gap.""" + +PCRASTER_MISSING = 1e31 +"""The PCRaster missing value. + +Written as ``Float32`` it does not round-trip exactly (``float(np.float32(1e31))`` +is ``9.999999...e30``), so the mask rejects everything at or above ``1e30`` +instead of comparing with this constant: no streamflow or runoff of a real +basin comes anywhere near that magnitude. +""" + +_MISSING_THRESHOLD = 1e30 + +INADMISSIBLE_OBJECTIVE = 1.0e30 +"""The objective of a candidate that is not run. + +The worst objective a run can produce is bounded by the efficiencies it can +reach: an NSE of -100, already an absurd simulation, gives +``1000 * (100 * 101) ** 2``, about ``1.0e11``. ``1.0e30`` is far above any such +value, so an inadmissible or failed candidate is always ranked behind every +candidate that was actually evaluated. It is finite so that it sorts, is +written to JSON and lands in the table of evaluations like any other objective. +""" + + +def valid_mask(values: np.ndarray) -> np.ndarray: + """Return the mask of the values that count as an observation. + + A value counts when it is finite, is not one of :data:`MISSING_VALUES` and + is below ``1e30``, which covers the PCRaster missing value in any of its + round-trips. + + :param values: Values of one station, in step order. + :type values: numpy.ndarray + + :return: A boolean array of the shape of ``values``. + :rtype: numpy.ndarray + """ + array = np.asarray(values, dtype=np.float64) + mask = np.isfinite(array) & (array < _MISSING_THRESHOLD) + for missing in MISSING_VALUES: + mask &= array != missing + return mask + + +def nash_sutcliffe(simulated: np.ndarray, observed: np.ndarray) -> float | None: + """Return the Nash-Sutcliffe efficiency of one station, or ``None``. + + The efficiency is computed over the pairs whose two members are valid, from + its definition:: + + NSE = 1 - sum((simulated - observed) ** 2) / sum((observed - mean(observed)) ** 2) + + Both sums run over the same pairs, and the mean is the mean of the observed + values of those pairs: a step dropped because the simulation has no value + for it is dropped from the mean as well, so that the numerator and the + denominator always describe the same sample. + + :param simulated: Simulated values, in step order. + :type simulated: numpy.ndarray + + :param observed: Observed values, aligned with ``simulated``. + :type observed: numpy.ndarray + + :return: The efficiency, or ``None`` when fewer than two pairs are valid or + when the valid observed values are all equal, which leaves the + denominator at zero and the efficiency undefined. + :rtype: float | None + + :raises ValueError: If the two series do not have the same length. + """ + simulated = np.asarray(simulated, dtype=np.float64) + observed = np.asarray(observed, dtype=np.float64) + if simulated.shape != observed.shape: + raise ValueError( + "The simulated and the observed series must have the same shape, got " + f"{simulated.shape} and {observed.shape}." + ) + + mask = valid_mask(simulated) & valid_mask(observed) + if int(np.count_nonzero(mask)) < 2: + return None + + simulated_valid = simulated[mask] + observed_valid = observed[mask] + # The zero-variance guard: a station whose observations are all equal has no + # variability for the efficiency to explain. Constancy is tested on the + # values themselves, since the sum of the squared deviations of a constant + # series is only approximately zero in floating point. + if np.all(observed_valid == observed_valid[0]): + return None + + residuals = float(np.sum((simulated_valid - observed_valid) ** 2)) + variability = float(np.sum((observed_valid - np.mean(observed_valid)) ** 2)) + return 1.0 - residuals / variability + + +def objective(nse: float) -> float: + """Return the objective value of an efficiency, ``1000 (100 (1 - NSE))^2``. + + The objective is zero for a perfect simulation and grows quadratically as + the efficiency falls, which is what the differential evolution minimizes. + + :param nse: The Nash-Sutcliffe efficiency of the candidate. + :type nse: float + + :return: The objective value. + :rtype: float + """ + return 1000.0 * (100.0 * (1.0 - nse)) ** 2 + + +@dataclass(frozen=True) +class Series: + """A time series table: the steps and the values of each station. + + :param steps: The time steps, one entry per row, as integers. + :type steps: numpy.ndarray + + :param stations: Station id to the values of that station, as ``float64`` + arrays aligned with :attr:`steps`. + :type stations: dict[str, numpy.ndarray] + """ + + steps: np.ndarray + stations: dict[str, np.ndarray] + + +def _as_value(token: str) -> float: + """Return ``token`` as a float, or ``NaN`` when it is empty or not a number. + + An entry a table leaves empty or writes in a form that is not a number is a + gap; it becomes ``NaN`` and the mask drops it, as it drops the gaps written + as ``-9999`` or as the PCRaster missing value. + """ + try: + return float(token) + except ValueError: + return float("nan") + + +def _as_step(token: str, path: Path, number: int) -> int: + """Return the time step written in the first column of a row. + + :raises ValueError: If the token is not a number, naming the file and the + row, since a table whose first column is a date or a label is a common + mistake and the bare conversion error does not say where it came from. + """ + try: + return int(float(token)) + except ValueError: + raise ValueError( + f"Row {number} of the time series file {path} does not start with a time " + f"step: {token!r} is not a number." + ) from None + + +def _build(ids: list[str], steps: list[int], rows: list[list[str]], path: Path) -> Series: + """Assemble a :class:`Series` from the parsed header, steps and value rows.""" + duplicated = sorted({station for station in ids if ids.count(station) > 1}) + if duplicated: + raise ValueError( + f"The time series file {path} has repeated station ids: {', '.join(duplicated)}." + ) + if not ids: + raise ValueError(f"The time series file {path} has no station column.") + if not rows: + raise ValueError(f"The time series file {path} has no data rows.") + + values = np.array( + [[_as_value(token) for token in row] for row in rows], + dtype=np.float64, + ) + return Series( + steps=np.asarray(steps, dtype=np.int64), + stations={station: values[:, column] for column, station in enumerate(ids)}, + ) + + +def _read_csv(path: Path) -> Series: + """Read the CSV layout of the model's time series tables. + + The separator is ``;``, the header is ``0;;...`` and every row + carries the time step in the first column and one value per station. + """ + with path.open(encoding="utf8", newline="") as table: + records = [ + row for row in csv.reader(table, delimiter=";") if any(cell.strip() for cell in row) + ] + if not records: + raise ValueError(f"The time series file {path} is empty.") + + header = [cell.strip() for cell in records[0]] + ids = header[1:] + steps: list[int] = [] + rows: list[list[str]] = [] + for number, record in enumerate(records[1:], 2): + cells = [cell.strip() for cell in record] + if len(cells) != len(header): + raise ValueError( + f"Row {number} of the time series file {path} has {len(cells)} column(s), " + f"the header gives {len(header)}." + ) + steps.append(_as_step(cells[0], path, number)) + rows.append(cells[1:]) + return _build(ids, steps, rows, path) + + +def _read_tss(path: Path) -> Series: + """Read a PCRaster time series (``.tss``) file. + + The layout is a title line, a line with the number of columns, one line per + column name starting with the time step column, and then one whitespace + separated row per time step. The model writes its own ``.tss`` files + without that header (see ``rubem.file._file_conversions.tss2csv``, which + reads every line of them as data); such a file is read as well, and its + stations are then numbered ``1``, ``2``, ... in column order, which is the + numbering of the sample locations of a usual dataset. + """ + lines = [line for line in path.read_text(encoding="utf8").splitlines() if line.strip()] + if not lines: + raise ValueError(f"The time series file {path} is empty.") + + ids: list[str] | None = None + data_start = 0 + if len(lines) > 1 and lines[1].strip().isdigit(): + columns = int(lines[1].strip()) + names = [line.strip() for line in lines[2 : 2 + columns]] + if len(names) == columns and columns >= 2: + ids = names[1:] + data_start = 2 + columns + + rows: list[list[str]] = [] + steps: list[int] = [] + for number, line in enumerate(lines[data_start:], data_start + 1): + cells = line.split() + if ids is None: + ids = [str(column) for column in range(1, len(cells))] + if len(cells) != len(ids) + 1: + raise ValueError( + f"Row {number} of the time series file {path} has {len(cells)} column(s), " + f"the column names give {len(ids) + 1}." + ) + steps.append(_as_step(cells[0], path, number)) + rows.append(cells[1:]) + return _build(ids or [], steps, rows, path) + + +def read_series(path: PathInput) -> Series: + """Read a time series table, in the model's CSV layout or as a PCRaster TSS. + + The layout is recognised from the content: a first line carrying the ``;`` + separator is the CSV table the model writes, anything else is read as a + PCRaster time series. Missing entries are kept as they are written; the + mask of :func:`valid_mask` is what decides whether a value counts. + + :param path: The file to read. + :type path: str | os.PathLike[str] + + :return: The parsed series. + :rtype: Series + + :raises FileNotFoundError: If the file does not exist. + :raises ValueError: If the file is empty, has no station column, repeats a + station id or has a row that does not match its header. + """ + file_path = as_path(path) + text = file_path.read_text(encoding="utf8") + first_line = next((line for line in text.splitlines() if line.strip()), "") + if ";" in first_line: + return _read_csv(file_path) + return _read_tss(file_path) + + +def evaluate_series( + simulated: Series, + observed: Series, + spinup_steps: int = 0, +) -> tuple[float, dict[str, float | None]]: + """Return the mean efficiency of a run and the efficiency of each station. + + The two series are aligned on the steps they share, after the steps up to + ``spinup_steps`` have been dropped from both, and on the station ids they + share: a station or a step present on one side only is ignored. The mean is + taken over the stations that have an efficiency, so that a station of zero + variance, reported as ``None``, does not drag the mean. + + A step repeated in a table is aligned on its first row, since that is the + occurrence :func:`numpy.intersect1d` reports for a repeated value. + + :param simulated: The series the run wrote. + :type simulated: Series + + :param observed: The observed series. + :type observed: Series + + :param spinup_steps: Number of initial steps to exclude, defaults to ``0``. + :type spinup_steps: int, optional + + :return: The mean efficiency and the efficiency of every shared station. + :rtype: tuple[float, dict[str, float | None]] + + :raises ValueError: If the two series share no station id, share no time + step, or if no shared station has an efficiency. + """ + shared_ids = sorted(set(simulated.stations) & set(observed.stations)) + if not shared_ids: + raise ValueError( + "The simulated and the observed series have no station in common: the " + f"simulation has {', '.join(sorted(simulated.stations)) or 'none'} and the " + f"observations have {', '.join(sorted(observed.stations)) or 'none'}." + ) + + simulated_steps = np.asarray(simulated.steps, dtype=np.int64) + observed_steps = np.asarray(observed.steps, dtype=np.int64) + simulated_kept = simulated_steps > spinup_steps + observed_kept = observed_steps > spinup_steps + shared_steps, simulated_rows, observed_rows = np.intersect1d( + simulated_steps[simulated_kept], + observed_steps[observed_kept], + return_indices=True, + ) + if shared_steps.size == 0: + raise ValueError( + "The simulated and the observed series have no time step in common after " + f"the {spinup_steps} spin-up step(s): the simulation has " + f"{simulated_steps[simulated_kept].tolist()} and the observations have " + f"{observed_steps[observed_kept].tolist()}." + ) + + simulated_index = np.flatnonzero(simulated_kept)[simulated_rows] + observed_index = np.flatnonzero(observed_kept)[observed_rows] + + station_nse: dict[str, float | None] = { + station: nash_sutcliffe( + simulated.stations[station][simulated_index], + observed.stations[station][observed_index], + ) + for station in shared_ids + } + + values = [nse for nse in station_nse.values() if nse is not None] + if not values: + raise ValueError( + "No station has a Nash-Sutcliffe efficiency: every shared station " + f"({', '.join(shared_ids)}) has fewer than two valid pairs or observations " + "without variance." + ) + return float(np.mean(values)), station_nse diff --git a/tests/unit/calibration/__init__.py b/tests/unit/calibration/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/unit/calibration/test_objective.py b/tests/unit/calibration/test_objective.py new file mode 100644 index 00000000..bc76f666 --- /dev/null +++ b/tests/unit/calibration/test_objective.py @@ -0,0 +1,396 @@ +"""The calibration objective: masks, Nash-Sutcliffe efficiency, tables and alignment.""" + +import numpy as np +import pytest + +from rubem.calibration.objective import ( + INADMISSIBLE_OBJECTIVE, + MISSING_VALUES, + PCRASTER_MISSING, + Series, + evaluate_series, + nash_sutcliffe, + objective, + read_series, + valid_mask, +) + +# The PCRaster missing value as it comes back from a Float32 raster or table: +# 1e31 has no exact Float32 representation, so what is read is slightly below it. +FLOAT32_MISSING = float(np.float32(1e31)) + +# Every gap marker, used wherever a test drops a pair through one of them. +GAPS = (float("nan"), -9999.0, 1e31, FLOAT32_MISSING) + + +def series(steps, **stations): + """A :class:`Series` built from plain lists, as the readers produce it.""" + return Series( + steps=np.asarray(steps, dtype=np.int64), + stations={ + station: np.asarray(values, dtype=np.float64) for station, values in stations.items() + }, + ) + + +class TestValidMask: + @pytest.mark.unit + def test_ordinary_values_are_valid(self): + assert valid_mask(np.array([0.0, -1.5, 3.25, 1e29])).tolist() == [True] * 4 + + @pytest.mark.unit + @pytest.mark.parametrize("gap", GAPS) + def test_every_gap_marker_is_masked(self, gap): + assert valid_mask(np.array([1.0, gap, 2.0])).tolist() == [True, False, True] + + @pytest.mark.unit + def test_the_documented_markers_are_the_ones_the_mask_rejects(self): + assert MISSING_VALUES == (-9999.0,) + assert PCRASTER_MISSING == 1e31 + assert not valid_mask(np.array([PCRASTER_MISSING])).any() + + @pytest.mark.unit + def test_infinities_are_masked(self): + assert valid_mask(np.array([np.inf, -np.inf])).tolist() == [False, False] + + @pytest.mark.unit + def test_the_threshold_itself_is_the_first_value_rejected(self): + # The mask keeps everything strictly below 1e30 and drops 1e30 itself, + # which is what makes the Float32 round-trip of 1e31 a gap as well. + below = float(np.nextafter(1e30, 0.0)) + + assert valid_mask(np.array([below, 1e30])).tolist() == [True, False] + assert below < 1e30 < FLOAT32_MISSING + + +class TestNashSutcliffe: + @pytest.mark.unit + def test_an_exact_simulation_gives_one(self): + # residuals 0 / variability 2 -> NSE = 1 - 0 = 1.0 + observed = np.array([1.0, 2.0, 3.0]) + + assert nash_sutcliffe(observed.copy(), observed) == 1.0 + + @pytest.mark.unit + def test_the_mean_of_the_observations_gives_zero(self): + # observed [1, 2, 3], mean 2; simulated [2, 2, 2] + # residuals (1 + 0 + 1) = 2, variability (1 + 0 + 1) = 2 -> NSE = 1 - 2/2 = 0.0 + assert nash_sutcliffe(np.array([2.0, 2.0, 2.0]), np.array([1.0, 2.0, 3.0])) == 0.0 + + @pytest.mark.unit + def test_a_reversed_simulation_gives_a_negative_efficiency(self): + # observed [1, 2, 3], mean 2; simulated [3, 2, 1] + # residuals (4 + 0 + 4) = 8, variability 2 -> NSE = 1 - 8/2 = -3.0 + assert nash_sutcliffe(np.array([3.0, 2.0, 1.0]), np.array([1.0, 2.0, 3.0])) == -3.0 + + @pytest.mark.unit + @pytest.mark.parametrize("gap", GAPS) + def test_a_gap_in_the_observations_drops_the_pair(self, gap): + # The fourth pair is dropped, the first three are exact -> NSE = 1.0 + simulated = np.array([1.0, 2.0, 3.0, 99.0]) + observed = np.array([1.0, 2.0, 3.0, gap]) + + assert nash_sutcliffe(simulated, observed) == 1.0 + + @pytest.mark.unit + @pytest.mark.parametrize("gap", GAPS) + def test_a_gap_in_the_simulation_drops_the_pair(self, gap): + # Dropping the fourth pair leaves observed [1, 2, 3] with mean 2 and + # simulated [3, 2, 1]: residuals 8, variability 2 -> NSE = -3.0 + simulated = np.array([3.0, 2.0, 1.0, gap]) + observed = np.array([1.0, 2.0, 3.0, 4.0]) + + assert nash_sutcliffe(simulated, observed) == -3.0 + + @pytest.mark.unit + def test_a_gap_is_dropped_from_the_mean_of_the_observations(self): + # The third pair is dropped by the simulation, so the mean is that of + # [1, 3] = 2 and not that of [1, 3, 8]: variability (1 + 1) = 2, + # residuals (1 + 1) = 2 -> NSE = 1 - 2/2 = 0.0 + simulated = np.array([2.0, 2.0, float("nan")]) + observed = np.array([1.0, 3.0, 8.0]) + + assert nash_sutcliffe(simulated, observed) == 0.0 + + @pytest.mark.unit + def test_a_single_valid_pair_has_no_efficiency(self): + assert nash_sutcliffe(np.array([1.0, 2.0]), np.array([1.0, float("nan")])) is None + + @pytest.mark.unit + def test_a_station_without_one_valid_pair_has_no_efficiency(self): + # No pair survives the mask at all: the count guard must answer before + # anything is read out of the empty selection. + assert nash_sutcliffe(np.array([1.0, -9999.0]), np.array([np.nan, 2.0])) is None + assert nash_sutcliffe(np.array([np.nan, np.nan]), np.array([np.nan, np.nan])) is None + + @pytest.mark.unit + def test_two_valid_pairs_are_enough(self): + # observed [1, 3], mean 2, variability 2; simulated [1, 3] -> NSE = 1.0 + simulated = np.array([1.0, 3.0, 5.0]) + observed = np.array([1.0, 3.0, -9999.0]) + + assert nash_sutcliffe(simulated, observed) == 1.0 + + @pytest.mark.unit + def test_observations_without_variance_have_no_efficiency(self): + # The sum of the squared deviations of [0.1, 0.1, 0.1] is not exactly + # zero in floating point, so constancy is what the guard must test. + assert nash_sutcliffe(np.array([0.1, 0.2, 0.3]), np.array([0.1, 0.1, 0.1])) is None + + @pytest.mark.unit + def test_only_the_valid_observations_decide_the_variance(self): + # The valid observations are [5, 5]: constant, so there is no efficiency + # although the series as written also carries a 7. + assert nash_sutcliffe(np.array([5.0, 5.0, np.nan]), np.array([5.0, 5.0, 7.0])) is None + + @pytest.mark.unit + def test_series_of_different_lengths_are_refused(self): + with pytest.raises(ValueError, match="same shape"): + nash_sutcliffe(np.array([1.0, 2.0]), np.array([1.0, 2.0, 3.0])) + + +class TestObjective: + @pytest.mark.unit + @pytest.mark.parametrize( + ("nse", "expected"), + [ + # 1000 * (100 * (1 - 1))^2 = 0 + (1.0, 0.0), + # 1000 * (100 * 1)^2 = 1000 * 10000 = 1e7 + (0.0, 1.0e7), + # 1000 * (100 * 2)^2 = 1000 * 40000 = 4e7 + (-1.0, 4.0e7), + ], + ) + def test_the_objective_of_a_known_efficiency(self, nse, expected): + assert objective(nse) == expected + + @pytest.mark.unit + def test_the_inadmissible_objective_is_far_above_any_run(self): + # An NSE of -100 is an absurd simulation and still gives ~1.02e11. + assert objective(-100.0) < INADMISSIBLE_OBJECTIVE / 1e10 + assert np.isfinite(INADMISSIBLE_OBJECTIVE) + + +class TestReadSeries: + @pytest.mark.unit + def test_a_handwritten_table_in_the_model_csv_layout(self, tmp_path): + path = tmp_path / "tss_arn.csv" + path.write_text("0;1;2\n1;1.5;2.5\n2;3.5;4.5\n3;5.5;6.5\n", encoding="utf8") + + table = read_series(path) + + assert table.steps.tolist() == [1, 2, 3] + assert sorted(table.stations) == ["1", "2"] + assert table.stations["1"].tolist() == [1.5, 3.5, 5.5] + assert table.stations["2"].tolist() == [2.5, 4.5, 6.5] + + @pytest.mark.unit + def test_gaps_of_a_csv_table_are_kept_as_written(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;A;B\n1;-9999;2.0\n2;1e31;\n", encoding="utf8") + + table = read_series(path) + + assert table.stations["A"].tolist() == [-9999.0, 1e31] + assert table.stations["B"][0] == 2.0 + assert np.isnan(table.stations["B"][1]), "an empty cell is a gap" + assert valid_mask(table.stations["A"]).tolist() == [False, False] + + @pytest.mark.unit + def test_the_ids_of_a_real_run_are_the_ids_of_the_sample_raster(self, tmp_path): + """The table the model writes must be readable exactly as it lands.""" + from tests.unit.core.test_core import run_model + + run_model(str(tmp_path)) + + table = read_series(tmp_path / "out" / "tss_arn.csv") + + # The synthetic dataset marks two sample locations, with the ids 1 and 2. + assert sorted(table.stations) == ["1", "2"] + assert table.steps.tolist() == [1, 2] + for values in table.stations.values(): + assert values.shape == (2,) + assert valid_mask(values).all() + + @pytest.mark.unit + def test_the_windows_line_endings_the_model_writes(self, tmp_path): + # rubem.file._file_conversions.tss2csv writes through csv.writer, whose + # line terminator is CRLF; the reader must not take the carriage return + # for part of the last value of a row. + path = tmp_path / "tss_arn.csv" + path.write_bytes(b"0;1;2\r\n1;1.5;2.5\r\n2;3.5;4.5\r\n") + + table = read_series(path) + + assert table.steps.tolist() == [1, 2] + assert table.stations["2"].tolist() == [2.5, 4.5] + + @pytest.mark.unit + def test_the_exponent_spellings_of_the_pcraster_missing_value(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;A\n1;1e+31\n2;1E31\n3;2.0\n", encoding="utf8") + + table = read_series(path) + + assert valid_mask(table.stations["A"]).tolist() == [False, False, True] + + @pytest.mark.unit + def test_a_handwritten_pcraster_time_series(self, tmp_path): + path = tmp_path / "observed.tss" + path.write_text( + "Observed streamflow\n3\ntimestep\n11\n22\n1 1.5 2.5\n2 3.5 4.5\n", + encoding="utf8", + ) + + table = read_series(path) + + assert table.steps.tolist() == [1, 2] + assert sorted(table.stations) == ["11", "22"] + assert table.stations["11"].tolist() == [1.5, 3.5] + assert table.stations["22"].tolist() == [2.5, 4.5] + + @pytest.mark.unit + def test_a_pcraster_time_series_without_the_header(self, tmp_path): + # The form the model itself writes: every line is a data row and the + # stations are numbered in column order. + path = tmp_path / "observed.tss" + path.write_text("1 1.5 2.5\n2 3.5 4.5\n", encoding="utf8") + + table = read_series(path) + + assert table.steps.tolist() == [1, 2] + assert sorted(table.stations) == ["1", "2"] + assert table.stations["2"].tolist() == [2.5, 4.5] + + @pytest.mark.unit + def test_a_row_that_does_not_start_with_a_time_step_names_the_file(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;A\n2000-01-01;1.0\n2000-01-02;2.0\n", encoding="utf8") + + with pytest.raises(ValueError, match="Row 2 .* does not start with a time step"): + read_series(path) + + @pytest.mark.unit + def test_a_repeated_station_id_is_refused(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;1;1\n1;1.0;2.0\n", encoding="utf8") + + with pytest.raises(ValueError, match="repeated station ids: 1"): + read_series(path) + + @pytest.mark.unit + def test_a_row_that_does_not_match_the_header_is_refused(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;1;2\n1;1.0;2.0\n2;3.0\n", encoding="utf8") + + with pytest.raises(ValueError, match="Row 3"): + read_series(path) + + @pytest.mark.unit + def test_an_empty_table_is_refused(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("", encoding="utf8") + + with pytest.raises(ValueError, match="empty"): + read_series(path) + + +class TestEvaluateSeries: + @pytest.mark.unit + def test_the_mean_over_the_stations(self): + # Station A: exact -> 1.0. Station B: observed [1, 2, 3] with the + # simulation reversed -> -3.0. Mean = (1.0 - 3.0) / 2 = -1.0 + simulated = series([1, 2, 3], A=[4.0, 5.0, 6.0], B=[3.0, 2.0, 1.0]) + observed = series([1, 2, 3], A=[4.0, 5.0, 6.0], B=[1.0, 2.0, 3.0]) + + mean, per_station = evaluate_series(simulated, observed) + + assert per_station == {"A": 1.0, "B": -3.0} + assert mean == -1.0 + + @pytest.mark.unit + def test_only_the_shared_steps_are_compared(self): + # The shared steps are 3 and 4: observed [1, 3] with mean 2 and + # variability 2, simulated [1, 3] -> NSE = 1.0 + simulated = series([1, 2, 3, 4], A=[99.0, 99.0, 1.0, 3.0]) + observed = series([3, 4, 5], A=[1.0, 3.0, 99.0]) + + mean, per_station = evaluate_series(simulated, observed) + + assert per_station == {"A": 1.0} + assert mean == 1.0 + + @pytest.mark.unit + def test_the_spin_up_steps_are_excluded(self): + # Without the spin-up the first step would make the simulation exact; + # with two spin-up steps only the steps 3 and 4 are compared: + # observed [1, 3], mean 2, variability 2; simulated [3, 1] + # -> residuals (4 + 4) = 8 -> NSE = 1 - 8/2 = -3.0 + simulated = series([1, 2, 3, 4], A=[0.0, 0.0, 3.0, 1.0]) + observed = series([1, 2, 3, 4], A=[0.0, 0.0, 1.0, 3.0]) + + mean, per_station = evaluate_series(simulated, observed, spinup_steps=2) + + assert per_station == {"A": -3.0} + assert mean == -3.0 + + @pytest.mark.unit + def test_a_station_present_on_one_side_only_is_ignored(self): + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[1.0, 2.0, 3.0]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0], C=[1.0, 2.0, 3.0]) + + mean, per_station = evaluate_series(simulated, observed) + + assert per_station == {"A": 1.0} + assert mean == 1.0 + + @pytest.mark.unit + def test_a_station_without_variance_is_reported_and_left_out_of_the_mean(self): + # Station A is exact -> 1.0; station B has constant observations -> None. + # The mean is that of the stations with a value: 1.0 + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[0.1, 0.2, 0.3]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[0.1, 0.1, 0.1]) + + mean, per_station = evaluate_series(simulated, observed) + + assert per_station == {"A": 1.0, "B": None} + assert mean == 1.0 + + @pytest.mark.unit + def test_no_shared_station_names_both_sides(self): + simulated = series([1, 2], A=[1.0, 2.0]) + observed = series([1, 2], B=[1.0, 2.0]) + + with pytest.raises(ValueError, match="no station in common") as error: + evaluate_series(simulated, observed) + + assert "A" in str(error.value) and "B" in str(error.value) + + @pytest.mark.unit + def test_no_shared_step_names_the_steps(self): + simulated = series([1, 2], A=[1.0, 2.0]) + observed = series([7, 8], A=[1.0, 2.0]) + + with pytest.raises(ValueError, match="no time step in common") as error: + evaluate_series(simulated, observed) + + assert "[1, 2]" in str(error.value) and "[7, 8]" in str(error.value) + + @pytest.mark.unit + def test_a_spin_up_that_eats_every_step_is_an_error(self): + simulated = series([1, 2], A=[1.0, 2.0]) + observed = series([1, 2], A=[1.0, 2.0]) + + with pytest.raises(ValueError, match="no time step in common"): + evaluate_series(simulated, observed, spinup_steps=2) + + @pytest.mark.unit + def test_no_station_with_an_efficiency_is_an_error(self): + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[1.0, 2.0, 3.0]) + observed = series([1, 2, 3], A=[5.0, 5.0, 5.0], B=[1.0, -9999.0, np.nan]) + + with pytest.raises(ValueError, match="No station has a Nash-Sutcliffe") as error: + evaluate_series(simulated, observed) + + assert "A, B" in str(error.value) From 0fb5b6d78c940166e2414dd2ea39cd23e03a05c7 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Wed, 16 Sep 2026 17:50:40 -0300 Subject: [PATCH 02/12] Add the differential evolution calibration runner parameters.py maps the eight free parameters to the decision vector, with w_3 derived from the other two weights, the bounds of the application settings and the linear constraint w_1 + w_2 <= 1. _worker.py is the module-level entry point a spawned worker runs: it builds the candidate's format 1.0 document with every raster series disabled and only the target time series enabled, runs the model in-process through the Python API, compares the table with the observed series and records one JSON file atomically; failures are recorded and ranked last. runner.py loads and validates the configuration once, drives SciPy's differential_evolution over a spawn ProcessPoolExecutor with one task per silenced worker, and writes evaluations.csv, result.json and the calibrated configuration. --- rubem/calibration/_worker.py | 281 +++++++++ rubem/calibration/parameters.py | 184 ++++++ rubem/calibration/runner.py | 678 ++++++++++++++++++++++ tests/unit/calibration/test_parameters.py | 227 ++++++++ tests/unit/calibration/test_runner.py | 580 ++++++++++++++++++ tests/unit/calibration/test_worker.py | 220 +++++++ 6 files changed, 2170 insertions(+) create mode 100644 rubem/calibration/_worker.py create mode 100644 rubem/calibration/parameters.py create mode 100644 rubem/calibration/runner.py create mode 100644 tests/unit/calibration/test_parameters.py create mode 100644 tests/unit/calibration/test_runner.py create mode 100644 tests/unit/calibration/test_worker.py diff --git a/rubem/calibration/_worker.py b/rubem/calibration/_worker.py new file mode 100644 index 00000000..edf32e89 --- /dev/null +++ b/rubem/calibration/_worker.py @@ -0,0 +1,281 @@ +"""Evaluation of one candidate, the entry point the calibration workers run. + +The module lives at the top level of the package and keeps its entry point, +:func:`evaluate`, at module level, because the process pool of the calibration +starts its workers with the ``spawn`` method: the child imports this module by +name to find the function it has to call, so the function must be picklable by +reference and everything it needs must travel as plain data in an +:class:`EvaluationContext`. + +One evaluation builds the configuration of the candidate from the document of +the calibrated configuration, runs the model in the current process, reads the +time series it wrote and compares it with the observed series. It never calls +:meth:`rubem.api.Model.run_isolated`: the worker is already the isolation layer, +and an isolated run inside it would nest one process pool in another. + +Nothing here raises: a candidate outside the admissible region is rejected +without a run, and any failure of the run or of the evaluation is recorded and +turned into :data:`rubem.calibration.objective.INADMISSIBLE_OBJECTIVE`, so that +one failing candidate does not end the calibration. +""" + +import copy +import json +import logging +import os +import shutil +import tempfile +import time +import uuid +from collections.abc import Sequence +from dataclasses import dataclass +from pathlib import Path + +import numpy as np + +from ..configuration.model_configuration_file_v1 import VARIABLE_IDS +from .objective import ( + INADMISSIBLE_OBJECTIVE, + Series, + evaluate_series, + objective, + read_series, +) +from .parameters import is_admissible, vector_to_parameters + +logger = logging.getLogger(__name__) + +_INADMISSIBLE_ERROR = "inadmissible" +"""The error recorded for a candidate that was rejected without being run.""" + + +@dataclass(frozen=True) +class EvaluationContext: + """Everything an evaluation needs, as data that crosses to a worker process. + + :param document: The calibrated configuration as a format 1.0 document, + with its paths already anchored. + :type document: dict + + :param base_dir: Directory the relative paths of the document are anchored + on, ``None`` when they are absolute already. + :type base_dir: str | None + + :param variable: Id of the output variable the efficiency is computed on + (``arn``, the accumulated total runoff, by default). + :type variable: str + + :param observed: The observed series, read once by the parent. + :type observed: rubem.calibration.objective.Series + + :param spinup_steps: Number of initial time steps excluded from the + efficiency. + :type spinup_steps: int + + :param temp_dir: Directory the per-evaluation output directories are + created in and removed from. + :type temp_dir: str + + :param evaluations_dir: Directory the JSON record of every evaluation is + written to. + :type evaluations_dir: str + + :param validate_input: Whether every evaluation revalidates the input files + and their content. ``False`` by default: the parent validated them once + before the calibration started, and they do not change during it. + :type validate_input: bool + """ + + document: dict + base_dir: str | None + variable: str + observed: Series + spinup_steps: int + temp_dir: str + evaluations_dir: str + validate_input: bool = False + + +def evaluate(vector: Sequence[float] | np.ndarray, context: EvaluationContext) -> float: + """Run one candidate and return the objective value of its simulation. + + A candidate outside the admissible region is rejected before any model run + and recorded with the error ``inadmissible``. Otherwise the model is run in + the current process, into a temporary output directory that is removed + again, and the efficiency of the resulting time series is turned into the + objective. Every evaluation, successful or not, leaves one JSON record in + the evaluations directory of the context. + + :param vector: The candidate, the eight free values in the order of + :data:`rubem.calibration.parameters.FREE_PARAMETERS`. + :type vector: collections.abc.Sequence[float] | numpy.ndarray + + :param context: What the evaluation runs against. + :type context: EvaluationContext + + :return: The objective value, or + :data:`rubem.calibration.objective.INADMISSIBLE_OBJECTIVE` when the + candidate was rejected or its evaluation failed. + :rtype: float + """ + started = time.perf_counter() + try: + parameters = vector_to_parameters(vector) + except ValueError as error: + _write_record(context, {}, None, {}, INADMISSIBLE_OBJECTIVE, 0.0, _describe(error)) + return INADMISSIBLE_OBJECTIVE + + if not is_admissible(vector): + logger.debug("Rejecting an inadmissible candidate without running the model.") + _write_record( + context, + parameters, + None, + {}, + INADMISSIBLE_OBJECTIVE, + time.perf_counter() - started, + _INADMISSIBLE_ERROR, + ) + return INADMISSIBLE_OBJECTIVE + + output_dir = tempfile.mkdtemp(prefix="evaluation-", dir=context.temp_dir) + try: + nse, station_nse = _run_and_evaluate(context, parameters, output_dir) + except Exception as error: + # One failing candidate must not end the calibration: it is recorded and + # ranked behind every candidate that was actually evaluated. + logger.warning("The evaluation of a candidate failed: %s", error) + _write_record( + context, + parameters, + None, + {}, + INADMISSIBLE_OBJECTIVE, + time.perf_counter() - started, + _describe(error), + ) + return INADMISSIBLE_OBJECTIVE + finally: + shutil.rmtree(output_dir, ignore_errors=True) + + value = objective(nse) + _write_record(context, parameters, nse, station_nse, value, time.perf_counter() - started, None) + return value + + +def _run_and_evaluate( + context: EvaluationContext, + parameters: dict[str, float], + output_dir: str, +) -> tuple[float, dict[str, float | None]]: + """Run the candidate and return its mean efficiency and the efficiency per station.""" + from ..api import Model + + document = _derived_document(context, parameters, output_dir) + result = Model.from_config( + document, validate_input=context.validate_input, base_dir=context.base_dir + ).run() + simulated = read_series(_time_series_table(result, context.variable)) + return evaluate_series(simulated, context.observed, context.spinup_steps) + + +def _time_series_table(result, variable: str) -> Path: + """Return the CSV table the run wrote for ``variable``. + + :raises RuntimeError: If the run reports no table for the variable. + """ + tables = result.time_series.get(variable, ()) + for path in tables: + if path.suffix.lower() == ".csv": + return path + raise RuntimeError( + f"The run wrote no CSV time series for the variable '{variable}'; it reports " + f"{', '.join(str(path) for path in tables) or 'no table'}." + ) + + +def _derived_document( + context: EvaluationContext, + parameters: dict[str, float], + output_dir: str, +) -> dict: + """Return the configuration document of one candidate. + + The calibrated document is copied and four things are changed: the + calibration parameters become the candidate's, the outputs go to the + temporary directory of the evaluation, every raster series is disabled (an + evaluation reads no raster, and writing them would dominate its cost) and + the time series are reduced to the calibrated variable, as CSV. The + aggregation the user configured is kept, since it decides which areas the + observed stations correspond to. + """ + document = copy.deepcopy(context.document) + document["model_calibration_parameters"] = { + "alpha": parameters["alpha"], + "b": parameters["beta"], + "w_1": parameters["w_1"], + "w_2": parameters["w_2"], + "w_3": parameters["w_3"], + "rcd": parameters["rcd"], + "f": parameters["f"], + "alpha_gw": parameters["alpha_gw"], + "x": parameters["x"], + } + output = document.setdefault("model_simulation_output", {}) + previous_rasters = output.get("raster_series", {}) + previous_samples = output.get("time_series_samples", {}) + output["dir_path"] = output_dir + output["raster_series"] = { + "no_data_value": previous_rasters.get("no_data_value", -9999), + "formats": [], + } + output["time_series_samples"] = { + **dict.fromkeys(VARIABLE_IDS, False), + context.variable: True, + "formats": ["CSV"], + "aggregation": previous_samples.get("aggregation", "point"), + } + return document + + +def _describe(error: BaseException) -> str: + """Return the recorded description of an exception: its type and its message.""" + return f"{type(error).__name__}: {error}" + + +def _write_record( + context: EvaluationContext, + parameters: dict[str, float], + nse: float | None, + station_nse: dict[str, float | None], + value: float, + elapsed_seconds: float, + error: str | None, +) -> None: + """Write the JSON record of one evaluation, atomically. + + Every worker writes into the same directory while the parent may be reading + it, so a record is written to a temporary file next to its destination and + renamed onto it: a reader either does not see the record yet or sees it + whole. The name is a fresh UUID, so two workers never collide. + """ + record = { + "id": uuid.uuid4().hex, + "pid": os.getpid(), + "parameters": parameters, + "nse": nse, + "station_nse": station_nse, + "objective": value, + "elapsed_seconds": elapsed_seconds, + "error": error, + } + directory = Path(context.evaluations_dir) + target = directory / f"{record['id']}.json" + handle, temporary = tempfile.mkstemp(prefix=f".{target.name}.", suffix=".tmp", dir=directory) + try: + with os.fdopen(handle, "w", encoding="utf-8") as file: + json.dump(record, file) + Path(temporary).replace(target) + except BaseException: + Path(temporary).unlink(missing_ok=True) + raise diff --git a/rubem/calibration/parameters.py b/rubem/calibration/parameters.py new file mode 100644 index 00000000..de14d6ac --- /dev/null +++ b/rubem/calibration/parameters.py @@ -0,0 +1,184 @@ +"""The decision vector of the calibration and its translation to model parameters. + +The differential evolution searches eight free parameters, ``alpha``, ``beta``, +``w_1``, ``w_2``, ``rcd``, ``f``, ``alpha_gw`` and ``x``; the slope factor +weight ``w_3`` is not searched but derived from the other two weights, since the +three weights of the weighted runoff coefficient must add up to 1. The search is +therefore eight-dimensional with the linear constraint ``w_1 + w_2 <= 1``, which +keeps the derived weight inside its own range. + +The bounds are the ones the application settings declare (``appsettings.json``, +section ``value_ranges.variables``), the same ranges +:class:`rubem.configuration.calibration_parameters.CalibrationParameters` +validates against, so a candidate the optimizer proposes inside the bounds is a +configuration the model accepts. + +Nothing here imports SciPy at module level; :func:`weights_constraint`, the only +function that needs it, imports it when it is called. +""" + +import logging +from collections.abc import Mapping, Sequence +from typing import TYPE_CHECKING + +import numpy as np + +from ..configuration._ranges import variable_range + +if TYPE_CHECKING: + from scipy.optimize import LinearConstraint + +logger = logging.getLogger(__name__) + +FREE_PARAMETERS = ("alpha", "beta", "w_1", "w_2", "rcd", "f", "alpha_gw", "x") +"""The parameters the optimizer searches, in the order of the decision vector.""" + +DERIVED_PARAMETER = "w_3" +"""The slope factor weight, derived from ``w_1`` and ``w_2`` instead of searched.""" + +CALIBRATION_PARAMETERS = ( + "alpha", + "beta", + "w_1", + "w_2", + "w_3", + "rcd", + "f", + "alpha_gw", + "x", +) +"""The nine calibration parameters of the model, in the order of the configuration.""" + +_ALIASES = {"b": "beta", "w1": "w_1", "w2": "w_2", "w3": "w_3"} +"""Spellings the legacy configuration file accepts for a calibration parameter.""" + + +def bounds() -> list[tuple[float, float]]: + """Return the ``(minimum, maximum)`` bound of every free parameter. + + The bounds come from the application settings, in the order of + :data:`FREE_PARAMETERS`, which is the order of the decision vector. + + :return: One ``(minimum, maximum)`` pair per free parameter. + :rtype: list[tuple[float, float]] + """ + return [variable_range(name) for name in FREE_PARAMETERS] + + +def vector_to_parameters(vector: Sequence[float] | np.ndarray) -> dict[str, float]: + """Return the nine calibration parameters a decision vector stands for. + + The eight free values are read in the order of :data:`FREE_PARAMETERS` and + the slope factor weight is derived as ``w_3 = 1 - w_1 - w_2``, so that the + three weights always add up to 1. + + :param vector: The eight free values, in the order of :data:`FREE_PARAMETERS`. + :type vector: collections.abc.Sequence[float] | numpy.ndarray + + :return: The nine parameters, keyed by the names of + :class:`rubem.configuration.calibration_parameters.CalibrationParameters`. + :rtype: dict[str, float] + + :raises ValueError: If the vector does not have eight entries. + """ + values = np.asarray(vector, dtype=np.float64).ravel() + if values.size != len(FREE_PARAMETERS): + raise ValueError( + f"The decision vector must have {len(FREE_PARAMETERS)} entries " + f"({', '.join(FREE_PARAMETERS)}), got {values.size}." + ) + parameters = {name: float(value) for name, value in zip(FREE_PARAMETERS, values, strict=True)} + # The two weights are added first and the sum is subtracted once. Taking + # them away one at a time, ``1 - w_1 - w_2``, rounds twice: for 20 of the + # 101 pairs of two-decimal weights that add up to 1, among them + # ``w_1 = 0.33, w_2 = 0.67``, it lands one ulp below zero, and + # :class:`rubem.configuration.calibration_parameters.CalibrationParameters` + # then refuses the run. With a single subtraction the weight is exactly + # zero at the boundary and never negative while ``w_1 + w_2 <= 1``. + parameters[DERIVED_PARAMETER] = 1.0 - (parameters["w_1"] + parameters["w_2"]) + return {name: parameters[name] for name in CALIBRATION_PARAMETERS} + + +def parameters_to_vector(parameters: Mapping[str, float]) -> np.ndarray: + """Return the decision vector of a set of calibration parameters. + + This is how the configuration under calibration becomes the ``x0`` of the + differential evolution. Both spellings the configuration files use are + accepted, the canonical ``beta``, ``w_1``, ``w_2`` and the legacy ``b``, + ``w1``, ``w2``; the derived ``w_3`` is ignored, since it is not searched. + The straightforward way to obtain the mapping is to dump the loaded + :class:`rubem.configuration.calibration_parameters.CalibrationParameters` + of a configuration, which has already normalised the spellings. + + :param parameters: The calibration parameters, by name. + :type parameters: collections.abc.Mapping[str, float] + + :return: The eight free values, in the order of :data:`FREE_PARAMETERS`. + :rtype: numpy.ndarray + + :raises KeyError: If a free parameter is missing from the mapping. + """ + normalised = {_ALIASES.get(name, name): value for name, value in parameters.items()} + missing = [name for name in FREE_PARAMETERS if name not in normalised] + if missing: + raise KeyError( + f"The calibration parameters are missing {', '.join(missing)}; " + f"the decision vector needs {', '.join(FREE_PARAMETERS)}." + ) + return np.asarray([float(normalised[name]) for name in FREE_PARAMETERS], dtype=np.float64) + + +def is_admissible(vector: Sequence[float] | np.ndarray) -> bool: + """Whether a decision vector describes a configuration the model accepts. + + A vector is admissible when every free value lies inside its bound, when + the two searched weights add up to at most 1 and when the derived ``w_3`` + lies inside its own bound. The derived weight is taken from + :func:`vector_to_parameters`, the very value the run would be given, so + that the guard can never disagree with the derivation and let through a + candidate + :class:`rubem.configuration.calibration_parameters.CalibrationParameters` + would then refuse. + + :param vector: The eight free values, in the order of :data:`FREE_PARAMETERS`. + :type vector: collections.abc.Sequence[float] | numpy.ndarray + + :return: ``True`` when the model can be run with this candidate. + :rtype: bool + """ + values = np.asarray(vector, dtype=np.float64).ravel() + if values.size != len(FREE_PARAMETERS): + return False + if not np.all(np.isfinite(values)): + return False + for value, (minimum, maximum) in zip(values, bounds(), strict=True): + if not minimum <= value <= maximum: + return False + + weights = float(values[FREE_PARAMETERS.index("w_1")] + values[FREE_PARAMETERS.index("w_2")]) + if weights > 1.0: + return False + derived = vector_to_parameters(values)[DERIVED_PARAMETER] + minimum, maximum = variable_range(DERIVED_PARAMETER) + return minimum <= derived <= maximum + + +def weights_constraint() -> "LinearConstraint": + """Return the linear constraint ``w_1 + w_2 <= 1`` on the decision vector. + + The differential evolution receives it so that it never spends a model run + on a candidate whose slope factor weight would come out negative: SciPy + evaluates the constraint before the objective and skips the members that + violate it. + + :return: The constraint, over the eight entries of the decision vector. + :rtype: scipy.optimize.LinearConstraint + + :raises ImportError: If SciPy is not installed. + """ + from scipy.optimize import LinearConstraint + + coefficients = np.zeros((1, len(FREE_PARAMETERS)), dtype=np.float64) + coefficients[0, FREE_PARAMETERS.index("w_1")] = 1.0 + coefficients[0, FREE_PARAMETERS.index("w_2")] = 1.0 + return LinearConstraint(coefficients, -np.inf, 1.0) diff --git a/rubem/calibration/runner.py b/rubem/calibration/runner.py new file mode 100644 index 00000000..4f99fb6f --- /dev/null +++ b/rubem/calibration/runner.py @@ -0,0 +1,678 @@ +"""The calibration itself: a differential evolution over the free parameters. + +:func:`calibrate` loads and validates the configuration once, reads the observed +series, and hands SciPy's ``differential_evolution`` a process pool whose workers +run one model evaluation each and exit. The pool is the isolation layer: PCRaster +keeps its clone and its raster memory process-wide, so a worker that is reused +would accumulate the state of every run it made; ``spawn`` with +``max_tasks_per_child=1`` gives every evaluation a fresh interpreter. + +The run directory receives three artifacts: ``evaluations.csv``, one row per +evaluation, ``result.json``, the best candidate and the settings that found it, +and ``-calibrated.json``, the calibrated configuration in the +format the input was written in. The rows of the table are ordered by the +candidate they evaluated and not by the moment they were written, so that two +calibrations of one configuration with one seed produce the same table. + +Number of evaluations + The population of a ``sobol`` initialization is ``max(5, popsize * 8)`` + rounded up to the next power of two, since a Sobol' sequence is balanced + only over a power-of-two sample: with the default ``popsize`` of 15 that is + 128 members, so a calibration costs one model run per member and per + generation, ``128 * (maxiter + 1)`` runs at most. The number is logged when + the calibration starts so that the budget can be checked before it runs. + Members whose weights violate ``w_1 + w_2 <= 1`` are skipped by SciPy + without a model run, so the real count is usually lower. + +SciPy is an optional dependency, imported inside :func:`calibrate`; importing +this module without it works and raises nothing. +""" + +import csv +import functools +import json +import logging +import math +import multiprocessing +import os +import sys +import tempfile +from concurrent.futures import ProcessPoolExecutor +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any + +import numpy as np + +from .._paths import PathInput, as_path +from ..configuration.model_configuration_file_v1 import ( + VARIABLE_IDS, + Aggregation, + ModelConfigurationFileV1, +) +from ._worker import _INADMISSIBLE_ERROR, EvaluationContext, evaluate +from .objective import read_series +from .parameters import ( + CALIBRATION_PARAMETERS, + FREE_PARAMETERS, + bounds, + parameters_to_vector, + vector_to_parameters, + weights_constraint, +) + +logger = logging.getLogger(__name__) + +EVALUATIONS_DIRNAME = "evaluations" +"""Name of the directory under the run directory that collects the JSON records.""" + +EVALUATIONS_CSV = "evaluations.csv" +RESULT_JSON = "result.json" + +EVALUATION_COLUMNS = ( + "id", + "pid", + *CALIBRATION_PARAMETERS, + "nse", + "objective", + "elapsed_seconds", + "error", +) +"""The columns of ``evaluations.csv``, one row per evaluation.""" + + +class CalibrationError(RuntimeError): + """A calibration that could not start, or that produced no usable evaluation.""" + + +@dataclass(frozen=True) +class CalibrationSettings: + """How the differential evolution is run. + + The differential evolution defaults follow the method the calibration of the + model was designed around; the remaining fields are the knobs the command + line exposes. + + :param variable: Id of the output variable compared with the observations. + ``arn``, the accumulated total runoff in m3/s, by default; ``rnf`` is + the total runoff per cell in mm. + :type variable: str + + :param spinup_steps: Number of initial time steps excluded from the + efficiency. + :type spinup_steps: int + + :param maxiter: Maximum number of generations. + :type maxiter: int + + :param popsize: Population multiplier; see the module docstring for the + number of members it produces. + :type popsize: int + + :param seed: Seed of the random generator of the search, for a reproducible + calibration. ``None`` leaves the generator unseeded. + :type seed: int | None + + :param workers: Number of worker processes. ``None``, the default, resolves + at call time to one less than the number of CPUs, at least one, which + leaves a core for the parent. + :type workers: int | None + + :param temp_dir: Directory the per-evaluation output directories are + created in. ``None`` uses the temporary directory of the system. + :type temp_dir: str | None + + :param init: Initial population: the name of a SciPy initialization + (``sobol`` by default) or an explicit ``(S, 8)`` array with ``S > 4``. + :type init: str | numpy.ndarray + + :param strategy: Differential evolution strategy. + :type strategy: str + + :param mutation: Mutation constant, a ``(min, max)`` pair for dithering. + :type mutation: tuple[float, float] + + :param recombination: Crossover probability. + :type recombination: float + + :param polish: Whether a local search refines the best candidate. ``False``: + the local search would spend model runs on gradients the objective, a + simulation read from a table, does not provide reliably. + :type polish: bool + """ + + variable: str = "arn" + spinup_steps: int = 0 + maxiter: int = 100 + popsize: int = 15 + seed: int | None = None + workers: int | None = None + temp_dir: str | None = None + init: str | np.ndarray = "sobol" + strategy: str = "best1exp" + mutation: tuple[float, float] = (0.5, 1.0) + recombination: float = 0.7 + polish: bool = False + + +@dataclass(frozen=True) +class CalibrationResult: + """What a finished calibration found and where it wrote it. + + :param best_parameters: The nine calibration parameters of the best + candidate, the derived ``w_3`` included. + :type best_parameters: dict[str, float] + + :param best_nse: The mean Nash-Sutcliffe efficiency of the best candidate, + ``None`` when no recorded evaluation reached one. + :type best_nse: float | None + + :param best_objective: The objective value of the best candidate. + :type best_objective: float + + :param evaluations: Number of objective evaluations the search spent. + :type evaluations: int + + :param generations: Number of generations the search ran. + :type generations: int + + :param success: Whether SciPy reports the search as successful. + :type success: bool + + :param message: The termination message of SciPy. + :type message: str + + :param run_dir: The run directory the artifacts were written to. + :type run_dir: pathlib.Path + + :param evaluations_csv: The table of every evaluation. + :type evaluations_csv: pathlib.Path + + :param result_json: The summary of the calibration. + :type result_json: pathlib.Path + + :param calibrated_config: The input configuration with the best parameters. + :type calibrated_config: pathlib.Path + """ + + best_parameters: dict[str, float] + best_nse: float | None + best_objective: float + evaluations: int + generations: int + success: bool + message: str + run_dir: Path + evaluations_csv: Path + result_json: Path + calibrated_config: Path + + +def _prepare_worker() -> None: + """Silence a worker process before it evaluates its candidate. + + The pool starts every worker as a fresh interpreter with no logging + configuration, so the warnings the loader emits for a configuration whose + inputs are not validated again would reach the terminal through the last + resort handler of :mod:`logging`, once per evaluation; the parent already + reported them when it validated the inputs. The PCRaster framework also + prints a newline when an interpreter exits, one per evaluation. A null + handler on the root logger and a standard output pointed at the null device + keep both out of the terminal: every failure of an evaluation is recorded in + its JSON record and summarized by the parent instead. + """ + logging.getLogger().addHandler(logging.NullHandler()) + sys.stdout = Path(os.devnull).open("w", encoding="utf-8") + + +def _report_failures(records: list[dict[str, Any]]) -> None: + """Log one summary of the evaluations that failed, if any did.""" + failures = [ + record["error"] + for record in records + if record.get("error") not in (None, _INADMISSIBLE_ERROR) + ] + if failures: + logger.warning( + "%d of %d evaluation(s) failed and were ranked behind every evaluated candidate; " + "the first failure was: %s. Every failure is recorded in the evaluations table.", + len(failures), + len(records), + failures[0], + ) + + +@dataclass +class _Progress: + """The per-generation callback: it only logs, it never stops the search.""" + + evaluations_dir: Path + generations: int = field(default=0) + + def __call__(self, intermediate_result) -> bool: + """Log the generation, the best objective so far and the records written. + + SciPy inspects the signature of the callback and passes the + intermediate result by keyword under exactly this name; returning + ``True`` would halt the search, so the callback always returns + ``False``. + """ + self.generations = int(getattr(intermediate_result, "nit", 0) or 0) + logger.info( + "Generation %d: best objective %.6g, %d evaluation(s) recorded.", + self.generations, + float(intermediate_result.fun), + len(_record_paths(self.evaluations_dir)), + ) + return False + + +def calibrate( + config_path: PathInput, + observed_path: PathInput, + run_dir: PathInput, + settings: CalibrationSettings = CalibrationSettings(), +) -> CalibrationResult: + """Calibrate a configuration against an observed series and write the artifacts. + + :param config_path: The configuration file to calibrate, legacy or format 1.0. + :type config_path: str | os.PathLike[str] + + :param observed_path: The observed series, in the CSV layout of the time + series tables of the model or as a PCRaster TSS file. + :type observed_path: str | os.PathLike[str] + + :param run_dir: Directory the artifacts of the calibration are written to. + It is created when it does not exist. + :type run_dir: str | os.PathLike[str] + + :param settings: How the search is run, defaults to :class:`CalibrationSettings`. + :type settings: CalibrationSettings, optional + + :return: The best candidate and the artifacts of the run. + :rtype: CalibrationResult + + :raises CalibrationError: If SciPy is not installed, if the configuration + cannot be calibrated as asked, or if every evaluation of the run failed. + :raises ImportError: If PCRaster or GDAL are not installed. + :raises FileNotFoundError: If the configuration or the observed series is + not there. + :raises ConfigurationError: If the configuration carries blocking problems. + """ + differential_evolution = _require_scipy() + + from ..api import Model + + if settings.variable not in VARIABLE_IDS: + raise CalibrationError( + f"'{settings.variable}' is not an output variable of the model; " + f"calibrate against one of {', '.join(VARIABLE_IDS)}." + ) + if settings.spinup_steps < 0: + raise CalibrationError( + f"The number of spin-up steps must not be negative, got {settings.spinup_steps}." + ) + if settings.workers is not None and settings.workers < 1: + raise CalibrationError( + f"The calibration needs at least one worker process, got {settings.workers}." + ) + + config_file = as_path(config_path).absolute() + configuration = Model.from_file(config_file, validate_input=True).configuration + file_v1 = _as_v1(configuration) + _check_sample_locations(configuration, file_v1) + + observed = read_series(observed_path) + x0 = parameters_to_vector(configuration.calibration_parameters.model_dump()) + + run_directory = as_path(run_dir).absolute() + evaluations_dir = run_directory / EVALUATIONS_DIRNAME + evaluations_dir.mkdir(parents=True, exist_ok=True) + previous = _record_paths(evaluations_dir) + if previous: + raise CalibrationError( + f"The run directory {run_directory} already holds {len(previous)} evaluation " + "record(s) of an earlier calibration; consolidating both runs into one table " + "would mix them. Choose an empty run directory." + ) + temp_dir = ( + as_path(settings.temp_dir).absolute() + if settings.temp_dir is not None + else Path(tempfile.gettempdir()) + ) + temp_dir.mkdir(parents=True, exist_ok=True) + + workers = ( + settings.workers if settings.workers is not None else max(1, (os.cpu_count() or 2) - 1) + ) + context = EvaluationContext( + document=file_v1.to_dict(), + base_dir=configuration.base_dir, + variable=settings.variable, + observed=observed, + spinup_steps=settings.spinup_steps, + temp_dir=str(temp_dir), + evaluations_dir=str(evaluations_dir), + validate_input=False, + ) + + planned_members = _population_size(settings) + logger.info( + "Calibrating %d free parameter(s) with %d population member(s) per generation and " + "at most %d generation(s): up to %d model run(s), on %d worker process(es).", + len(FREE_PARAMETERS), + planned_members, + settings.maxiter, + planned_members * (settings.maxiter + 1), + workers, + ) + + progress = _Progress(evaluations_dir) + executor = ProcessPoolExecutor( + max_workers=workers, + mp_context=multiprocessing.get_context("spawn"), + max_tasks_per_child=1, + initializer=_prepare_worker, + ) + try: + optimum = differential_evolution( + functools.partial(evaluate, context=context), + bounds(), + strategy=settings.strategy, + maxiter=settings.maxiter, + popsize=settings.popsize, + mutation=settings.mutation, + recombination=settings.recombination, + rng=settings.seed, + polish=settings.polish, + init=settings.init, + x0=x0, + updating="deferred", + workers=executor.map, + constraints=weights_constraint(), + callback=progress, + ) + finally: + executor.shutdown(wait=True) + + records = _read_records(evaluations_dir) + evaluations_csv = run_directory / EVALUATIONS_CSV + _write_evaluations_csv(evaluations_csv, records) + + _report_failures(records) + best = _best_record(records) + if best is None: + raise CalibrationError( + "Every evaluation of the calibration failed; the first failure was: " + f"{_first_error(records)}. The table of the evaluations is {evaluations_csv}." + ) + + # The size of the population SciPy actually built, not the size this module + # predicted for the budget: a named initialization resizes the population + # (``sobol`` rounds it up to a power of two) and the summary must describe + # the search that ran. + population_size = int(np.shape(optimum.population)[0]) + if population_size != planned_members: + logger.info( + "The search used %d population member(s), %d were announced.", + population_size, + planned_members, + ) + + best_parameters = vector_to_parameters(optimum.x) + best_objective = float(optimum.fun) + result_json = run_directory / RESULT_JSON + calibrated_config = run_directory / f"{config_file.stem}-calibrated.json" + _write_result_json( + result_json, + best_parameters=best_parameters, + best_nse=best["nse"], + best_objective=best_objective, + optimum=optimum, + settings=settings, + workers=workers, + temp_dir=temp_dir, + population_size=population_size, + ) + _write_calibrated_config(calibrated_config, configuration, file_v1, best_parameters) + + logger.info( + "Calibration finished after %d evaluation(s) in %d generation(s): best objective %.6g.", + int(optimum.nfev), + int(optimum.nit), + best_objective, + ) + return CalibrationResult( + best_parameters=best_parameters, + best_nse=best["nse"], + best_objective=best_objective, + evaluations=int(optimum.nfev), + generations=int(optimum.nit), + success=bool(optimum.success), + message=str(optimum.message), + run_dir=run_directory, + evaluations_csv=evaluations_csv, + result_json=result_json, + calibrated_config=calibrated_config, + ) + + +def _require_scipy(): + """Return ``scipy.optimize.differential_evolution``, or explain how to install it. + + :raises CalibrationError: If SciPy is not installed. + """ + from .. import _deps + + missing = _deps.missing_calibration_deps() + if missing: + raise CalibrationError(_deps.calibration_deps_message(missing)) + from scipy.optimize import differential_evolution + + return differential_evolution + + +def _as_v1(configuration) -> ModelConfigurationFileV1: + """Return the configuration as a format 1.0 model, converting a legacy one.""" + if configuration.file_v1 is not None: + return configuration.file_v1 + return ModelConfigurationFileV1.from_legacy(configuration.file) + + +def _check_sample_locations(configuration, file_v1: ModelConfigurationFileV1) -> None: + """Refuse a configuration whose time series would have nowhere to be sampled. + + :raises CalibrationError: If the raster the configured aggregation needs is + not part of the configuration. + """ + aggregation = file_v1.model_simulation_output.time_series_samples.aggregation + if aggregation is Aggregation.ZONES: + if not configuration.raster_files.zones: + raise CalibrationError( + "The configuration aggregates its time series over zones but names no " + "zones raster; the calibration has no station to compare." + ) + return + if not configuration.raster_files.sample_locations: + raise CalibrationError( + "The configuration names no sample locations raster; the calibration has no " + "station to compare the observed series with." + ) + + +def _population_size(settings: CalibrationSettings) -> int: + """Return the number of population members the search will use. + + The arithmetic mirrors SciPy's: an explicit population is used as it is, a + named initialization gives ``max(5, popsize * 8)`` members, and ``sobol`` + rounds that up to the next power of two. + """ + if not isinstance(settings.init, str): + return int(np.shape(settings.init)[0]) + members = max(5, settings.popsize * len(FREE_PARAMETERS)) + if settings.init == "sobol": + return 1 << (members - 1).bit_length() + return members + + +def _record_paths(evaluations_dir: Path) -> list[Path]: + """Return the JSON records written so far, without the files still being staged.""" + return sorted(path for path in evaluations_dir.glob("*.json") if not path.name.startswith(".")) + + +def _read_records(evaluations_dir: Path) -> list[dict[str, Any]]: + """Read every JSON record of the run, in the order of :func:`_record_order`. + + A record that cannot be read is reported and skipped: an unreadable record + of one evaluation must not hide the result of all the others. + """ + records = [] + for path in _record_paths(evaluations_dir): + try: + records.append(json.loads(path.read_text(encoding="utf-8"))) + except (OSError, ValueError) as error: + logger.warning("Skipping the unreadable evaluation record %s: %s", path, error) + records.sort(key=_record_order) + return records + + +def _record_order(record: dict[str, Any]) -> tuple[tuple[float, ...], str]: + """Return the sort key of a record: its parameters, then its identifier. + + The records are not ordered by the moment they were written: the workers run + in parallel, so that order changes from run to run. Ordering them by the + candidate they evaluated makes ``evaluations.csv`` a function of the seed + alone, so two calibrations of the same configuration with the same seed + produce the same table and can be compared line by line. The identifier only + separates two evaluations of the very same candidate. + """ + parameters = record.get("parameters") or {} + values = tuple(float(parameters.get(name, math.inf)) for name in CALIBRATION_PARAMETERS) + return values, str(record.get("id", "")) + + +def _write_evaluations_csv(path: Path, records: list[dict[str, Any]]) -> None: + """Write one row per evaluation, in the columns of :data:`EVALUATION_COLUMNS`.""" + with path.open("w", encoding="utf-8", newline="") as table: + writer = csv.writer(table) + writer.writerow(EVALUATION_COLUMNS) + for record in records: + parameters = record.get("parameters") or {} + writer.writerow( + [ + record.get("id", ""), + record.get("pid", ""), + *(_cell(parameters.get(name)) for name in CALIBRATION_PARAMETERS), + _cell(record.get("nse")), + _cell(record.get("objective")), + _cell(record.get("elapsed_seconds")), + record.get("error") or "", + ] + ) + + +def _cell(value) -> str: + """Return a numeric cell with full precision, or an empty cell for a missing value.""" + return "" if value is None else repr(float(value)) + + +def _best_record(records: list[dict[str, Any]]) -> dict[str, Any] | None: + """Return the record of the lowest objective among those that were evaluated.""" + evaluated = [record for record in records if not record.get("error")] + if not evaluated: + return None + return min(evaluated, key=lambda record: float(record["objective"])) + + +def _first_error(records: list[dict[str, Any]]) -> str: + """Return the error of the first record that carries one.""" + return next((str(record["error"]) for record in records if record.get("error")), "unknown") + + +def _write_result_json( + path: Path, + *, + best_parameters: dict[str, float], + best_nse: float | None, + best_objective: float, + optimum, + settings: CalibrationSettings, + workers: int, + temp_dir: Path, + population_size: int, +) -> None: + """Write the summary of the calibration.""" + document = { + "best_parameters": best_parameters, + "best_nse": best_nse, + "best_objective": best_objective, + "nfev": int(optimum.nfev), + "nit": int(optimum.nit), + "success": bool(optimum.success), + "message": str(optimum.message), + "population_size": population_size, + "settings": _settings_document(settings, workers, temp_dir), + } + path.write_text(json.dumps(document, indent=2) + "\n", encoding="utf-8") + + +def _settings_document( + settings: CalibrationSettings, workers: int, temp_dir: Path +) -> dict[str, Any]: + """Return the settings as JSON-able data, resolved the way the run used them. + + The number of workers and the temporary directory are the values the run + actually used, not the ``None`` the caller may have left them at, so that + the summary describes a calibration that can be repeated. + """ + init = settings.init if isinstance(settings.init, str) else np.asarray(settings.init).tolist() + return { + "variable": settings.variable, + "spinup_steps": settings.spinup_steps, + "maxiter": settings.maxiter, + "popsize": settings.popsize, + "seed": settings.seed, + "workers": workers, + "temp_dir": str(temp_dir), + "init": init, + "strategy": settings.strategy, + "mutation": list(settings.mutation), + "recombination": settings.recombination, + "polish": settings.polish, + } + + +def _write_calibrated_config( + path: Path, + configuration, + file_v1: ModelConfigurationFileV1, + best_parameters: dict[str, float], +) -> None: + """Write the calibrated configuration in the format the input was written in. + + The document goes through the configuration models, so what is written is a + configuration the loader accepts. Its paths are the ones the loader resolved, + which are absolute even when the input file wrote them relative to its own + directory. + """ + calibrated = file_v1.model_copy( + update={ + "model_calibration_parameters": type(file_v1.model_calibration_parameters)( + alpha=best_parameters["alpha"], + b=best_parameters["beta"], + w_1=best_parameters["w_1"], + w_2=best_parameters["w_2"], + w_3=best_parameters["w_3"], + rcd=best_parameters["rcd"], + f=best_parameters["f"], + alpha_gw=best_parameters["alpha_gw"], + x=best_parameters["x"], + ) + } + ) + document = ( + calibrated.to_dict() + if configuration.file_v1 is not None + else calibrated.to_legacy().to_dict() + ) + path.write_text(json.dumps(document, indent=2) + "\n", encoding="utf-8") diff --git a/tests/unit/calibration/test_parameters.py b/tests/unit/calibration/test_parameters.py new file mode 100644 index 00000000..90d0db87 --- /dev/null +++ b/tests/unit/calibration/test_parameters.py @@ -0,0 +1,227 @@ +"""The decision vector of the calibration: bounds, conversions and admissibility.""" + +import numpy as np +import pytest + +from rubem.calibration.parameters import ( + CALIBRATION_PARAMETERS, + DERIVED_PARAMETER, + FREE_PARAMETERS, + bounds, + is_admissible, + parameters_to_vector, + vector_to_parameters, + weights_constraint, +) +from rubem.configuration.app_settings import AppSettings +from rubem.configuration.calibration_parameters import CalibrationParameters +from rubem.configuration.model_configuration import ModelConfiguration +from tests.helpers.synthetic import write_synthetic_dataset + +# The ranges of rubem/appsettings.json, written out so that a change of the +# settings file, or of the order of the decision vector, has to be made here +# too instead of passing unnoticed through a comparison with the same source. +SETTINGS_BOUNDS = { + "alpha": (0.01, 10.0), + "beta": (0.01, 1.0), + "w_1": (0.0, 1.0), + "w_2": (0.0, 1.0), + "w_3": (0.0, 1.0), + "rcd": (1.0, 10.0), + "f": (0.01, 1.0), + "alpha_gw": (0.01, 1.0), + "x": (0.0, 1.0), +} + + +def admissible_vector(): + """A decision vector well inside every bound, with weights that add up below 1.""" + return np.array([4.5, 0.5, 0.333, 0.333, 5.0, 0.5, 0.5, 0.5], dtype=np.float64) + + +class TestBounds: + @pytest.mark.unit + def test_the_bounds_are_the_ranges_of_the_application_settings(self): + ranges = AppSettings.default().value_ranges.variables + + assert bounds() == [(ranges[name].min, ranges[name].max) for name in FREE_PARAMETERS] + + @pytest.mark.unit + def test_the_bounds_are_the_ones_written_in_the_settings_file(self): + assert bounds() == [SETTINGS_BOUNDS[name] for name in FREE_PARAMETERS] + assert [name for name in SETTINGS_BOUNDS if name != DERIVED_PARAMETER] == list( + FREE_PARAMETERS + ) + + @pytest.mark.unit + def test_the_free_parameters_do_not_include_the_derived_weight(self): + assert DERIVED_PARAMETER not in FREE_PARAMETERS + assert set(CALIBRATION_PARAMETERS) == set(FREE_PARAMETERS) | {DERIVED_PARAMETER} + assert len(bounds()) == len(FREE_PARAMETERS) == 8 + + +class TestConversions: + @pytest.mark.unit + def test_a_vector_becomes_the_nine_parameters_with_the_derived_weight(self): + parameters = vector_to_parameters(admissible_vector()) + + assert list(parameters) == list(CALIBRATION_PARAMETERS) + assert parameters["beta"] == 0.5 + assert parameters["w_3"] == 1.0 - (0.333 + 0.333) + assert parameters["w_1"] + parameters["w_2"] + parameters["w_3"] == pytest.approx(1.0) + + @pytest.mark.unit + def test_a_vector_of_the_wrong_size_is_rejected(self): + with pytest.raises(ValueError, match="must have 8 entries"): + vector_to_parameters([1.0, 2.0, 3.0]) + + @pytest.mark.unit + def test_the_configuration_parameters_round_trip_through_the_vector(self, tmp_path): + config = write_synthetic_dataset(str(tmp_path)) + configuration = ModelConfiguration(config, validate_input=False) + parameters = configuration.calibration_parameters.model_dump() + + vector = parameters_to_vector(parameters) + + assert vector_to_parameters(vector) == pytest.approx(parameters, abs=1e-12) + + @pytest.mark.unit + def test_the_legacy_spellings_of_the_parameters_are_accepted(self): + legacy = { + "alpha": 4.5, + "b": 0.5, + "w1": 0.333, + "w2": 0.333, + "w3": 0.334, + "rcd": 5.0, + "f": 0.5, + "alpha_gw": 0.5, + "x": 0.5, + } + + assert parameters_to_vector(legacy) == pytest.approx(admissible_vector()) + + @pytest.mark.unit + def test_a_missing_parameter_names_what_the_vector_needs(self): + with pytest.raises(KeyError, match="alpha_gw"): + parameters_to_vector({"alpha": 4.5, "beta": 0.5}) + + +class TestAdmissibility: + @pytest.mark.unit + def test_a_vector_inside_the_bounds_is_admissible(self): + assert is_admissible(admissible_vector()) + + @pytest.mark.unit + def test_a_vector_on_the_bounds_is_admissible(self): + vector = np.array([minimum for minimum, _ in bounds()], dtype=np.float64) + + assert is_admissible(vector) + + @pytest.mark.unit + @pytest.mark.parametrize("index", range(8)) + def test_a_value_above_its_bound_is_rejected(self, index): + vector = admissible_vector() + vector[index] = bounds()[index][1] + 1.0 + + assert not is_admissible(vector) + + @pytest.mark.unit + @pytest.mark.parametrize("index", range(8)) + def test_a_value_below_its_bound_is_rejected(self, index): + vector = admissible_vector() + vector[index] = bounds()[index][0] - 1.0 + + assert not is_admissible(vector) + + @pytest.mark.unit + def test_weights_that_add_up_above_one_are_rejected(self): + vector = admissible_vector() + vector[FREE_PARAMETERS.index("w_1")] = 0.6 + vector[FREE_PARAMETERS.index("w_2")] = 0.6 + + # Both weights are inside their own [0, 1] bound; what makes the + # candidate inadmissible is their sum, and the negative w_3 it derives. + assert vector_to_parameters(vector)["w_3"] < 0.0 + assert not is_admissible(vector) + + @pytest.mark.unit + def test_weights_that_add_up_to_exactly_one_are_admissible(self): + vector = admissible_vector() + vector[FREE_PARAMETERS.index("w_1")] = 0.4 + vector[FREE_PARAMETERS.index("w_2")] = 0.6 + + assert vector_to_parameters(vector)["w_3"] == 0.0 + assert is_admissible(vector) + + @pytest.mark.unit + @pytest.mark.parametrize("hundredths", range(101)) + def test_every_pair_of_weights_that_adds_up_to_one_reaches_the_model(self, hundredths): + """What the guard admits, the configuration model must accept. + + ``1 - w_1 - w_2`` rounds twice and lands one ulp below zero for 20 of + these 101 pairs, ``w_1 = 0.33, w_2 = 0.67`` among them, which is a + weight :class:`CalibrationParameters` refuses; deriving the weight with + a single subtraction gives exactly zero for all of them. + """ + vector = admissible_vector() + vector[FREE_PARAMETERS.index("w_1")] = hundredths / 100.0 + vector[FREE_PARAMETERS.index("w_2")] = (100 - hundredths) / 100.0 + parameters = vector_to_parameters(vector) + + assert parameters["w_3"] == 0.0 + assert is_admissible(vector) + # No exception: the candidate the guard admits is one the model runs. + CalibrationParameters(**parameters) + + @pytest.mark.unit + def test_the_guard_and_the_derivation_cannot_disagree(self): + """Every admissible candidate of a grid is one the model accepts.""" + refused = [] + for first in range(0, 101, 7): + for second in range(0, 101, 11): + vector = admissible_vector() + vector[FREE_PARAMETERS.index("w_1")] = first / 100.0 + vector[FREE_PARAMETERS.index("w_2")] = second / 100.0 + if not is_admissible(vector): + continue + try: + CalibrationParameters(**vector_to_parameters(vector)) + except ValueError: + refused.append((first / 100.0, second / 100.0)) + + assert not refused, f"admitted but refused by the model: {refused}" + + @pytest.mark.unit + def test_a_vector_of_the_wrong_size_or_with_a_nan_is_rejected(self): + vector = admissible_vector() + vector[0] = np.nan + + assert not is_admissible(vector) + assert not is_admissible(admissible_vector()[:-1]) + + +class TestWeightsConstraint: + """The only part of the module that needs the optional SciPy extra.""" + + @pytest.fixture(autouse=True) + def _scipy(self): + pytest.importorskip("scipy.optimize") + + @pytest.mark.unit + def test_the_constraint_is_the_sum_of_the_two_searched_weights(self): + constraint = weights_constraint() + + expected = np.zeros((1, len(FREE_PARAMETERS))) + expected[0, FREE_PARAMETERS.index("w_1")] = 1.0 + expected[0, FREE_PARAMETERS.index("w_2")] = 1.0 + assert np.array_equal(np.asarray(constraint.A), expected) + assert np.all(np.asarray(constraint.lb) == -np.inf) + assert np.all(np.asarray(constraint.ub) == 1.0) + + @pytest.mark.unit + def test_the_constraint_evaluates_to_the_sum_of_the_weights(self): + constraint = weights_constraint() + vector = admissible_vector() + + assert np.asarray(constraint.A) @ vector == pytest.approx([0.666]) diff --git a/tests/unit/calibration/test_runner.py b/tests/unit/calibration/test_runner.py new file mode 100644 index 00000000..d1fe1b53 --- /dev/null +++ b/tests/unit/calibration/test_runner.py @@ -0,0 +1,580 @@ +"""The calibration runner: a small deterministic differential evolution and its artifacts.""" + +import csv +import inspect +import json +import logging +import multiprocessing +import os +import shutil +import subprocess +import sys +from concurrent.futures import ProcessPoolExecutor +from pathlib import Path + +import numpy as np +import pytest + +from rubem import _deps +from rubem.api import Model +from rubem.calibration import runner +from rubem.calibration.parameters import ( + CALIBRATION_PARAMETERS, + FREE_PARAMETERS, + bounds, + parameters_to_vector, +) +from rubem.calibration.runner import ( + EVALUATION_COLUMNS, + CalibrationError, + CalibrationSettings, + _population_size, + _Progress, + calibrate, +) +from rubem.configuration.model_configuration import ModelConfiguration +from rubem.configuration.model_configuration_file import ModelConfigurationFile +from rubem.configuration.model_configuration_file_v1 import ModelConfigurationFileV1 +from tests.helpers.config import REPO_ROOT +from tests.helpers.synthetic import write_synthetic_dataset + +# SciPy is the optional ``rubem[calibration]`` extra, and an environment +# without it (the frozen one of the byte-exact job, for instance) still +# collects this file; the whole module is skipped there instead of failing to +# import. +scipy_optimize = pytest.importorskip("scipy.optimize") + +TIMESTEPS = 3 + +# Five admissible members, the second one far from the configuration under +# calibration. SciPy replaces the first member with x0, so the configuration +# itself is always part of the initial population. +INIT = np.array( + [ + [1.0, 0.20, 0.20, 0.20, 2.0, 0.20, 0.20, 0.20], + [9.5, 0.90, 0.10, 0.80, 9.0, 0.90, 0.90, 0.90], + [2.0, 0.30, 0.40, 0.40, 3.0, 0.30, 0.30, 0.30], + [6.0, 0.70, 0.25, 0.25, 7.0, 0.70, 0.70, 0.70], + [3.0, 0.40, 0.50, 0.30, 4.0, 0.40, 0.40, 0.40], + ] +) + + +class Dataset: + """The synthetic dataset, its configuration file and the series a plain run wrote.""" + + def __init__(self, tmp_path): + self.config = write_synthetic_dataset(str(tmp_path), timesteps=TIMESTEPS) + self.config_file = tmp_path / "config.json" + self.config_file.write_text(json.dumps(self.config), encoding="utf8") + written = Model.from_config(self.config).run().time_series["arn"][0] + self.observed = tmp_path / "observed.csv" + shutil.copyfile(written, self.observed) + self.run_dir = tmp_path / "calibration" + self.temp_dir = tmp_path / "temp" + self.parameters = ModelConfiguration( + self.config, validate_input=False + ).calibration_parameters.model_dump() + + def settings(self, **overrides): + """The deterministic small search of the tests.""" + defaults = { + "variable": "arn", + "maxiter": 1, + "popsize": 5, + "seed": 1, + "workers": 2, + "temp_dir": str(self.temp_dir), + "init": INIT, + "polish": False, + } + return CalibrationSettings(**{**defaults, **overrides}) + + def calibrate(self, observed=None, **overrides): + return calibrate( + self.config_file, + observed or self.observed, + self.run_dir, + self.settings(**overrides), + ) + + +def read_evaluations(path): + with path.open(encoding="utf-8", newline="") as table: + return list(csv.DictReader(table)) + + +@pytest.fixture +def dataset(tmp_path): + return Dataset(tmp_path) + + +@pytest.fixture(scope="class") +def calibrated(tmp_path_factory): + """One calibration, shared by the assertions that inspect its artifacts.""" + data = Dataset(tmp_path_factory.mktemp("calibration")) + return data, data.calibrate() + + +class TestCalibration: + @pytest.mark.unit + def test_the_configuration_that_produced_the_observations_is_the_optimum(self, calibrated): + data, result = calibrated + + # The observed series is what the configuration itself wrote, so the + # global optimum of the search is the configuration, with NSE 1. + assert result.best_objective == pytest.approx(0.0, abs=1e-9) + assert result.best_nse == pytest.approx(1.0, abs=1e-12) + assert result.best_parameters == pytest.approx(data.parameters, abs=1e-12) + assert set(result.best_parameters) == set(CALIBRATION_PARAMETERS) + + @pytest.mark.unit + def test_every_evaluation_is_a_row_of_the_table(self, calibrated): + data, result = calibrated + rows = read_evaluations(result.evaluations_csv) + + assert list(rows[0]) == list(EVALUATION_COLUMNS) + assert len(rows) == result.evaluations + assert result.best_objective == pytest.approx(min(float(row["objective"]) for row in rows)) + assert {row["id"] for row in rows} == {row["id"] for row in rows if row["id"]} + assert len(rows) >= len(INIT) + # The linear constraint keeps the candidates whose weights do not add up + # away from the workers, so no evaluation is ever rejected as + # inadmissible and none of them failed. + assert [row["error"] for row in rows if row["error"]] == [] + + @pytest.mark.unit + def test_the_evaluations_ran_in_other_processes_and_none_lingers(self, calibrated): + _, result = calibrated + rows = read_evaluations(result.evaluations_csv) + + pids = [int(row["pid"]) for row in rows] + + assert pids, "at least one evaluation was recorded" + assert os.getpid() not in pids + # One task per child: PCRaster keeps its clone and its rasters for the + # life of the process, so no two evaluations may share one. + assert len(set(pids)) == len(pids) + assert multiprocessing.active_children() == [] + + @pytest.mark.unit + def test_the_summary_carries_the_settings_and_the_population_size(self, calibrated): + data, result = calibrated + summary = json.loads(result.result_json.read_text(encoding="utf-8")) + + assert summary["nfev"] == result.evaluations + assert summary["nit"] == result.generations + assert summary["population_size"] == len(INIT) + assert summary["best_objective"] == pytest.approx(result.best_objective) + assert summary["best_nse"] == pytest.approx(result.best_nse) + assert summary["best_parameters"] == pytest.approx(result.best_parameters) + assert summary["settings"]["seed"] == 1 + assert summary["settings"]["workers"] == 2 + assert summary["settings"]["variable"] == "arn" + assert summary["settings"]["maxiter"] == 1 + assert summary["settings"]["polish"] is False + assert summary["settings"]["temp_dir"] == str(data.temp_dir) + assert summary["settings"]["init"] == INIT.tolist() + + @pytest.mark.unit + def test_the_calibrated_configuration_loads_with_the_best_parameters(self, calibrated): + data, result = calibrated + + assert result.calibrated_config.name == "config-calibrated.json" + configuration = ModelConfiguration(result.calibrated_config, validate_input=False) + assert configuration.calibration_parameters.model_dump() == pytest.approx( + result.best_parameters, abs=1e-12 + ) + document = json.loads(result.calibrated_config.read_text(encoding="utf-8")) + assert "CALIBRATION" in document, "a legacy configuration stays legacy" + + @pytest.mark.unit + def test_the_best_row_of_the_table_is_what_was_written_back(self, calibrated): + _, result = calibrated + rows = read_evaluations(result.evaluations_csv) + best = min(rows, key=lambda row: float(row["objective"])) + parameters = {name: float(best[name]) for name in CALIBRATION_PARAMETERS} + + # The optimum SciPy reports and the best record of the run are one and + # the same evaluation, down to the last bit, and it is that evaluation + # the calibrated configuration carries. + assert float(best["objective"]) == result.best_objective + assert float(best["nse"]) == result.best_nse + assert parameters == result.best_parameters + configuration = ModelConfiguration(result.calibrated_config, validate_input=False) + assert configuration.calibration_parameters.model_dump() == pytest.approx( + parameters, abs=1e-12 + ) + + @pytest.mark.unit + def test_the_temporary_directory_is_left_empty(self, calibrated): + data, _ = calibrated + + assert list(data.temp_dir.iterdir()) == [] + + @pytest.mark.unit + def test_the_run_directory_holds_the_three_artifacts(self, calibrated): + _, result = calibrated + + assert result.run_dir.is_dir() + assert result.evaluations_csv.is_file() + assert result.result_json.is_file() + assert result.calibrated_config.is_file() + assert result.message + assert result.generations == 1 + + +class _RecordingExecutor(ProcessPoolExecutor): + """A process pool that remembers how it was built and when it was closed.""" + + built: list["_RecordingExecutor"] = [] + + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.built_with = kwargs + self.shutdowns = 0 + _RecordingExecutor.built.append(self) + + def shutdown(self, *args, **kwargs): + self.shutdowns += 1 + super().shutdown(*args, **kwargs) + + +@pytest.fixture +def broken_search(monkeypatch): + """Record what the search is given, and let it fail without running a model.""" + call = {} + + def explode(function, search_bounds, **kwargs): + call["function"] = function + call["bounds"] = search_bounds + call.update(kwargs) + raise RuntimeError("the search broke") + + _RecordingExecutor.built.clear() + monkeypatch.setattr(scipy_optimize, "differential_evolution", explode) + monkeypatch.setattr(runner, "ProcessPoolExecutor", _RecordingExecutor) + return call + + +class TestProcessTopology: + @pytest.mark.unit + def test_the_pool_is_a_spawned_one_and_is_closed_when_the_search_fails( + self, dataset, broken_search + ): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate() + + (executor,) = _RecordingExecutor.built + assert executor.built_with["max_workers"] == 2 + assert executor.built_with["max_tasks_per_child"] == 1 + assert executor.built_with["mp_context"].get_start_method() == "spawn" + # The pool is closed even though the search raised, and the map of that + # very pool is what the search was told to evaluate with. + assert executor.shutdowns == 1 + assert broken_search["workers"].__self__ is executor + assert multiprocessing.active_children() == [] + + @pytest.mark.unit + def test_the_search_is_given_the_documented_arguments(self, dataset, broken_search): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate() + + assert broken_search["bounds"] == bounds() + assert broken_search["strategy"] == "best1exp" + assert broken_search["maxiter"] == 1 + assert broken_search["popsize"] == 5 + assert broken_search["mutation"] == (0.5, 1.0) + assert broken_search["recombination"] == 0.7 + assert broken_search["updating"] == "deferred" + assert broken_search["polish"] is False + assert broken_search["rng"] == 1 + assert broken_search["init"] is INIT + # The configuration under calibration is the first member of the + # population: SciPy puts x0 in place of the first row of init. + assert broken_search["x0"] == pytest.approx(parameters_to_vector(dataset.parameters)) + + @pytest.mark.unit + def test_the_search_is_constrained_to_weights_that_add_up_to_at_most_one( + self, dataset, broken_search + ): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate() + + constraint = broken_search["constraints"] + coefficients = np.atleast_2d(constraint.A) + expected = [ + [1.0 if name in ("w_1", "w_2") else 0.0 for name in FREE_PARAMETERS], + ] + + assert coefficients.tolist() == expected + assert np.atleast_1d(constraint.ub).tolist() == [1.0] + assert np.all(np.isneginf(np.atleast_1d(constraint.lb))) + + +class TestProgress: + @pytest.mark.unit + def test_the_callback_takes_the_argument_scipy_passes_by_keyword(self): + # SciPy reads the signature of the callback: only a callback whose one + # parameter is called ``intermediate_result`` is given the result of the + # generation, any other is called with the old ``(x, convergence)`` pair. + assert set(inspect.signature(_Progress(Path())).parameters) == {"intermediate_result"} + + @pytest.mark.unit + def test_the_callback_logs_the_generation_and_lets_the_search_go_on( + self, tmp_path, caplog, capfd + ): + evaluations = tmp_path / "evaluations" + evaluations.mkdir() + (evaluations / "one.json").write_text("{}", encoding="utf8") + progress = _Progress(evaluations) + + with caplog.at_level(logging.INFO, logger="rubem.calibration.runner"): + stop = progress(scipy_optimize.OptimizeResult(fun=12.5, nit=3)) + + assert stop is False, "the callback never halts the search" + assert progress.generations == 3 + assert "Generation 3" in caplog.text + assert "12.5" in caplog.text + assert "1 evaluation" in caplog.text + assert capfd.readouterr().out == "", "the library logs, it does not print" + + +class TestReproducibility: + @pytest.mark.unit + def test_two_runs_of_one_seed_write_the_same_table(self, dataset, capfd): + first = calibrate( + dataset.config_file, + dataset.observed, + dataset.run_dir / "first", + dataset.settings(), + ) + second = calibrate( + dataset.config_file, + dataset.observed, + dataset.run_dir / "second", + dataset.settings(), + ) + + def candidates(result): + return [ + tuple(row[name] for name in CALIBRATION_PARAMETERS) + for row in read_evaluations(result.evaluations_csv) + ] + + # Same seed, same initial population, same candidates, in the same order: + # the table of a calibration is a function of its settings alone. Only + # the identifiers and the process ids differ between the two runs. + assert candidates(first) == candidates(second) + assert first.best_parameters == second.best_parameters + assert first.best_objective == second.best_objective + assert first.evaluations == second.evaluations + # The library logs, the command line prints, and the workers are silent: + # the newline the PCRaster framework prints when a worker exits goes to + # the null device, and the warnings of an unvalidated load, which the + # parent already reported, meet a null handler instead of the last + # resort handler of logging. Captured at file-descriptor level, so the + # children are included. + out, err = capfd.readouterr() + assert out == "" + assert "WARNING" not in err and "No handlers could be found" not in err + + +class TestFailures: + @pytest.mark.unit + def test_a_run_whose_evaluations_all_fail_is_reported_with_its_first_error( + self, dataset, tmp_path + ): + foreign = tmp_path / "foreign.csv" + foreign.write_text("0;A;B\n1;1.0;2.0\n2;2.0;3.0\n3;3.0;4.0\n", encoding="utf8") + + with pytest.raises(CalibrationError, match="no station in common"): + dataset.calibrate(observed=foreign, maxiter=1) + + rows = read_evaluations(dataset.run_dir / "evaluations.csv") + assert rows, "the table of the evaluations is written even when every one failed" + assert all(row["error"] for row in rows) + assert all(row["nse"] == "" for row in rows) + + @pytest.mark.unit + def test_a_variable_the_model_does_not_write_is_refused(self, dataset): + with pytest.raises(CalibrationError, match="not an output variable"): + dataset.calibrate(variable="flow") + + @pytest.mark.unit + def test_a_run_directory_of_an_earlier_calibration_is_refused(self, dataset): + evaluations = dataset.run_dir / "evaluations" + evaluations.mkdir(parents=True) + (evaluations / "deadbeef.json").write_text("{}", encoding="utf8") + + with pytest.raises(CalibrationError, match="already holds 1 evaluation"): + dataset.calibrate() + + @pytest.mark.unit + def test_a_negative_spin_up_is_refused(self, dataset): + with pytest.raises(CalibrationError, match="must not be negative"): + dataset.calibrate(spinup_steps=-1) + + @pytest.mark.unit + def test_a_search_without_a_worker_is_refused(self, dataset): + with pytest.raises(CalibrationError, match="at least one worker"): + dataset.calibrate(workers=0) + + @pytest.mark.unit + def test_without_scipy_the_optional_extra_is_named(self, dataset, monkeypatch): + monkeypatch.setattr(_deps, "missing_calibration_deps", lambda: ["scipy"]) + + with pytest.raises(CalibrationError, match=r"rubem\[calibration\]"): + dataset.calibrate() + + assert not dataset.run_dir.exists(), "nothing is written when the search cannot start" + + +class TestSettings: + @pytest.mark.unit + def test_the_default_decision_vector_is_the_configuration_under_calibration(self, dataset): + vector = parameters_to_vector(dataset.parameters) + + assert vector == pytest.approx([4.5, 0.5, 0.333, 0.333, 5.0, 0.5, 0.5, 0.5]) + + @pytest.mark.unit + def test_the_defaults_are_the_documented_ones(self): + settings = CalibrationSettings() + + assert (settings.variable, settings.spinup_steps) == ("arn", 0) + assert (settings.maxiter, settings.popsize) == (100, 15) + assert settings.init == "sobol" + assert settings.strategy == "best1exp" + assert settings.mutation == (0.5, 1.0) + assert settings.recombination == 0.7 + assert settings.polish is False + assert settings.workers is None + + +@pytest.fixture(scope="class") +def calibrated_v1(tmp_path_factory): + """One calibration of the same dataset written as a format 1.0 configuration.""" + tmp_path = tmp_path_factory.mktemp("calibration-v1") + data = Dataset(tmp_path) + legacy = ModelConfigurationFile.model_validate(data.config) + data.config_file = tmp_path / "basin.json" + data.config_file.write_text( + json.dumps(ModelConfigurationFileV1.from_legacy(legacy).to_dict()), encoding="utf8" + ) + return data, data.calibrate() + + +class TestFormatV1Configuration: + @pytest.mark.unit + def test_a_format_1_0_configuration_is_calibrated_and_written_back_in_its_format( + self, calibrated_v1 + ): + data, result = calibrated_v1 + + assert result.best_objective == pytest.approx(0.0, abs=1e-9) + assert result.best_parameters == pytest.approx(data.parameters, abs=1e-12) + assert result.calibrated_config.name == "basin-calibrated.json" + document = json.loads(result.calibrated_config.read_text(encoding="utf-8")) + assert document["version"] == "1.0" + assert document["model_calibration_parameters"]["b"] == pytest.approx( + result.best_parameters["beta"], abs=1e-12 + ) + configuration = ModelConfiguration(result.calibrated_config, validate_input=False) + assert configuration.calibration_parameters.model_dump() == pytest.approx( + result.best_parameters, abs=1e-12 + ) + + +class TestEvaluationBudget: + @pytest.mark.unit + def test_a_sobol_population_is_rounded_up_to_a_power_of_two(self): + # The documented arithmetic: max(5, popsize * 8) members, rounded up to + # the next power of two, because a Sobol' sequence is balanced only over + # a power-of-two sample. + assert _population_size(CalibrationSettings()) == 128 + assert _population_size(CalibrationSettings(popsize=1)) == 8 + assert _population_size(CalibrationSettings(popsize=16)) == 128 + + @pytest.mark.unit + def test_another_initialization_uses_the_members_it_asks_for(self): + assert _population_size(CalibrationSettings(init="latinhypercube", popsize=3)) == 24 + assert _population_size(CalibrationSettings(init="random", popsize=15)) == 120 + assert _population_size(CalibrationSettings(init=INIT)) == len(INIT) + + @pytest.mark.unit + @pytest.mark.parametrize("popsize", [1, 5, 15, 16]) + @pytest.mark.parametrize("init", ["sobol", "latinhypercube", "random"]) + def test_the_announced_population_is_the_one_scipy_builds(self, popsize, init): + # The announced number is what the command line reports as the budget, + # so it has to be SciPy's own arithmetic and not an approximation of it. + optimum = scipy_optimize.differential_evolution( + lambda vector: float(np.sum(vector)), + bounds(), + init=init, + popsize=popsize, + maxiter=0, + polish=False, + rng=0, + ) + + assert len(optimum.population) == _population_size( + CalibrationSettings(popsize=popsize, init=init) + ) + + @pytest.mark.unit + def test_the_summary_reports_the_population_the_search_really_built(self, dataset, monkeypatch): + # A named initialization resizes the population, so the summary must + # report the population of the result and not the number this module + # announced for the budget. + def stub(function, search_bounds, **kwargs): + population = np.asarray(kwargs["init"], dtype=np.float64) + energies = np.asarray([function(member) for member in population[:2]]) + enlarged = np.repeat(population, 2, axis=0) + return scipy_optimize.OptimizeResult( + x=population[0], + fun=float(energies[0]), + nfev=len(energies), + nit=1, + success=True, + message="stub", + population=enlarged, + population_energies=np.full(len(enlarged), np.inf), + ) + + monkeypatch.setattr(scipy_optimize, "differential_evolution", stub) + + result = dataset.calibrate() + summary = json.loads(result.result_json.read_text(encoding="utf-8")) + + assert summary["population_size"] == 2 * len(INIT) + assert _population_size(dataset.settings()) == len(INIT) + + +class TestOptionalDependency: + @pytest.mark.unit + def test_importing_the_package_does_not_import_scipy(self, tmp_path): + # An installation without the rubem[calibration] extra must still be + # able to import the package, so SciPy is imported inside the functions + # that need it. A fresh interpreter is the only honest probe: the test + # session has SciPy loaded already. + probe = tmp_path / "probe.py" + probe.write_text( + "import sys\n" + "import rubem.calibration.runner\n" + "import rubem.calibration._worker\n" + "import rubem.calibration.parameters\n" + "print('scipy' in sys.modules)\n", + encoding="utf8", + ) + + environment = {**os.environ, "PYTHONPATH": REPO_ROOT} + completed = subprocess.run( + [sys.executable, str(probe)], + capture_output=True, + text=True, + check=True, + cwd=REPO_ROOT, + env=environment, + ) + + assert completed.stdout.strip() == "False" diff --git a/tests/unit/calibration/test_worker.py b/tests/unit/calibration/test_worker.py new file mode 100644 index 00000000..301d1011 --- /dev/null +++ b/tests/unit/calibration/test_worker.py @@ -0,0 +1,220 @@ +"""One evaluation of a candidate: the model run, the record and the failure modes.""" + +import json +import os +import shutil +from pathlib import Path + +import pytest + +from rubem.api import Model +from rubem.calibration._worker import EvaluationContext, _derived_document, evaluate +from rubem.calibration.objective import INADMISSIBLE_OBJECTIVE, read_series +from rubem.calibration.parameters import ( + FREE_PARAMETERS, + parameters_to_vector, + vector_to_parameters, +) +from rubem.configuration.model_configuration import ModelConfiguration +from rubem.configuration.model_configuration_file_v1 import ( + VARIABLE_IDS, + ModelConfigurationFileV1, +) +from tests.helpers.synthetic import write_synthetic_dataset + +TIMESTEPS = 3 + + +class Dataset: + """The synthetic dataset, its configuration and the series a plain run wrote.""" + + def __init__(self, tmp_path): + self.config = write_synthetic_dataset(str(tmp_path), timesteps=TIMESTEPS) + self.observed_path = Model.from_config(self.config).run().time_series["arn"][0] + self.configuration = ModelConfiguration(self.config, validate_input=False) + self.document = ModelConfigurationFileV1.from_legacy(self.configuration.file).to_dict() + self.vector = parameters_to_vector(self.configuration.calibration_parameters.model_dump()) + self.temp_dir = tmp_path / "temp" + self.evaluations_dir = tmp_path / "evaluations" + self.temp_dir.mkdir() + self.evaluations_dir.mkdir() + + def context(self, observed_path=None, variable="arn", spinup_steps=0): + """An evaluation context reading ``observed_path``, the plain run by default.""" + return EvaluationContext( + document=self.document, + base_dir=self.configuration.base_dir, + variable=variable, + observed=read_series(observed_path or self.observed_path), + spinup_steps=spinup_steps, + temp_dir=str(self.temp_dir), + evaluations_dir=str(self.evaluations_dir), + validate_input=False, + ) + + def records(self): + """Every record written so far, without the files still being staged.""" + return [ + json.loads(path.read_text(encoding="utf-8")) + for path in sorted(self.evaluations_dir.glob("*.json")) + if not path.name.startswith(".") + ] + + def staged(self): + """The staging files left behind by an interrupted write, if any.""" + return [path.name for path in self.evaluations_dir.iterdir() if ".tmp" in path.name] + + +@pytest.fixture +def dataset(tmp_path): + return Dataset(tmp_path) + + +class TestSuccessfulEvaluation: + @pytest.mark.unit + def test_the_configuration_of_the_observations_scores_a_perfect_objective(self, dataset): + # The observed series is what the configuration itself produced, so its + # own parameters reproduce it exactly: NSE 1, objective 0. + value = evaluate(dataset.vector, dataset.context()) + + assert value == pytest.approx(0.0, abs=1e-9) + (record,) = dataset.records() + assert record["error"] is None + assert record["nse"] == pytest.approx(1.0, abs=1e-12) + assert record["objective"] == pytest.approx(0.0, abs=1e-9) + assert record["station_nse"] == {"1": pytest.approx(1.0), "2": pytest.approx(1.0)} + assert record["parameters"]["w_3"] == pytest.approx(0.334, abs=1e-12) + assert record["pid"] == os.getpid() + assert record["elapsed_seconds"] >= 0.0 + + @pytest.mark.unit + def test_the_temporary_output_directory_is_removed_after_the_run(self, dataset): + evaluate(dataset.vector, dataset.context()) + + assert list(dataset.temp_dir.iterdir()) == [] + + @pytest.mark.unit + def test_the_record_is_installed_by_a_rename_and_leaves_no_staging_file(self, dataset): + evaluate(dataset.vector, dataset.context()) + + assert dataset.staged() == [] + assert len(dataset.records()) == 1 + + @pytest.mark.unit + def test_a_different_candidate_is_evaluated_and_scores_worse(self, dataset): + vector = dataset.vector.copy() + vector[FREE_PARAMETERS.index("alpha_gw")] = 0.01 + + value = evaluate(vector, dataset.context()) + + assert 0.0 < value < INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert record["error"] is None + assert record["nse"] < 1.0 + + +class TestDerivedConfiguration: + @pytest.mark.unit + def test_it_disables_every_raster_series_and_keeps_one_time_series(self, dataset): + document = _derived_document( + dataset.context(), vector_to_parameters(dataset.vector), "/output" + ) + output = document["model_simulation_output"] + rasters = output["raster_series"] + samples = output["time_series_samples"] + + assert output["dir_path"] == "/output" + assert rasters["formats"] == [] + assert [name for name in VARIABLE_IDS if rasters.get(name)] == [] + assert [name for name in VARIABLE_IDS if samples.get(name)] == ["arn"] + assert samples["formats"] == ["CSV"] + assert samples["aggregation"] == "point" + # What the worker hands to the model has to be a configuration again. + ModelConfigurationFileV1.model_validate(document) + + @pytest.mark.unit + def test_the_run_writes_nothing_but_the_table_of_the_calibrated_variable( + self, dataset, monkeypatch + ): + # The output directory is listed on its way out, so that the listing is + # of the directory the evaluation really wrote and the removal still + # happens exactly as it does without the test. + listed = [] + removal = shutil.rmtree + + def spy(path, **kwargs): + listed.extend(sorted(entry.name for entry in Path(path).rglob("*"))) + removal(path, **kwargs) + + monkeypatch.setattr(shutil, "rmtree", spy) + + evaluate(dataset.vector, dataset.context()) + + assert listed == ["metadata.json", "tss_arn.csv"] + + +class TestRejectedCandidate: + @pytest.mark.unit + def test_an_inadmissible_candidate_is_recorded_without_a_model_run(self, dataset): + vector = dataset.vector.copy() + vector[FREE_PARAMETERS.index("w_1")] = 0.6 + vector[FREE_PARAMETERS.index("w_2")] = 0.6 + + value = evaluate(vector, dataset.context()) + + assert value == INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert record["error"] == "inadmissible" + assert record["nse"] is None + assert record["objective"] == INADMISSIBLE_OBJECTIVE + assert list(dataset.temp_dir.iterdir()) == [], "no output directory was created" + + +class TestRecord: + @pytest.mark.unit + def test_the_record_is_staged_in_the_directory_it_is_renamed_into(self, dataset, monkeypatch): + # A rename is atomic only inside one filesystem, so the staging file has + # to be a sibling of the record and not a file in the system temporary + # directory, which may be a different mount. + renames = [] + rename = Path.replace + + def spy(self, target): + renames.append((Path(self), Path(target))) + return rename(self, target) + + monkeypatch.setattr(Path, "replace", spy) + + evaluate(dataset.vector, dataset.context()) + + installed = [pair for pair in renames if pair[1].parent == dataset.evaluations_dir] + assert len(installed) == 1, "the record was installed by exactly one rename" + temporary, target = installed[0] + assert temporary.parent == dataset.evaluations_dir + assert temporary.name.startswith(".") + assert target.name.endswith(".json") + + +class TestFailedEvaluation: + @pytest.mark.unit + def test_observations_of_another_station_are_recorded_as_the_error(self, dataset, tmp_path): + foreign = tmp_path / "foreign.csv" + foreign.write_text("0;A;B\n1;1.0;2.0\n2;2.0;3.0\n3;3.0;4.0\n", encoding="utf8") + + value = evaluate(dataset.vector, dataset.context(observed_path=foreign)) + + assert value == INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert "no station in common" in record["error"] + assert record["error"].startswith("ValueError: ") + assert record["nse"] is None + assert list(dataset.temp_dir.iterdir()) == [], "the output directory was removed anyway" + + @pytest.mark.unit + def test_a_variable_the_run_does_not_write_is_recorded_as_the_error(self, dataset): + value = evaluate(dataset.vector, dataset.context(variable="nope")) + + assert value == INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert record["error"] + assert list(dataset.temp_dir.iterdir()) == [] From 0da37fea5b202d89b32b005a9d2c08a219288110 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Wed, 16 Sep 2026 17:50:40 -0300 Subject: [PATCH 03/12] Add the rubem calibrate command The command takes the configuration, the observed series and the run directory, exposes the search settings as options, prints the documented progress lines and exits without a traceback on configuration and calibration errors. rubem._deps gains the optional calibration dependency guard, with the guidance to install rubem[calibration]. --- rubem/_deps.py | 50 ++++ rubem/cli.py | 114 +++++++- tests/unit/test_cli_calibrate.py | 479 +++++++++++++++++++++++++++++++ tests/unit/test_deps.py | 51 ++++ 4 files changed, 691 insertions(+), 3 deletions(-) create mode 100644 tests/unit/test_cli_calibrate.py diff --git a/rubem/_deps.py b/rubem/_deps.py index a27c0548..d6e846e1 100644 --- a/rubem/_deps.py +++ b/rubem/_deps.py @@ -69,3 +69,53 @@ def require_preprocessing_deps() -> None: "This preprocessing tool needs the optional dependencies " f"{', '.join(missing)}. Install them with 'pip install \"rubem[preprocessing]\"'." ) + + +_CALIBRATION_DEPENDENCIES = ("scipy",) + + +def missing_calibration_deps() -> list[str]: + """Return the optional calibration dependencies that are not importable. + + As in :func:`missing_runtime_deps`, a dependency whose lookup itself fails + counts as missing: the calibration would not get past it either. + + :return: The missing dependency names, in the order they are checked. + """ + missing = [] + for name in _CALIBRATION_DEPENDENCIES: + try: + found = importlib.util.find_spec(name) is not None + except (ImportError, ValueError): + found = False + if not found: + missing.append(name) + return missing + + +def calibration_deps_message(missing: list[str] | None = None) -> str: + """Return the installation guidance for the missing calibration dependencies. + + The command line raises it as ``SystemExit`` and the runner as + ``CalibrationError``, so both front ends give the same instructions. + + :param missing: The missing dependency names, see + :func:`missing_calibration_deps`. Defaults to ``None``, which names + every optional calibration dependency. + """ + names = list(_CALIBRATION_DEPENDENCIES) if missing is None else missing + plural = "dependencies" if len(names) > 1 else "dependency" + return ( + f"The calibration command needs the optional {plural} {', '.join(names)}. " + "Install it with 'pip install \"rubem[calibration]\"'." + ) + + +def require_calibration_deps() -> None: + """Fail fast with guidance when SciPy is not importable. + + :raises SystemExit: If any optional calibration dependency is missing. + """ + missing = missing_calibration_deps() + if missing: + raise SystemExit(calibration_deps_message(missing)) diff --git a/rubem/cli.py b/rubem/cli.py index fb9eb5fe..a1e0044e 100644 --- a/rubem/cli.py +++ b/rubem/cli.py @@ -1,6 +1,7 @@ """The ``rubem`` command line. -``rubem run -c [-s]`` runs a simulation; ``rubem config schema`` +``rubem run -c [-s]`` runs a simulation; ``rubem calibrate`` +fits the calibration parameters to an observed series; ``rubem config schema`` prints the JSON Schema of the configuration file. The former ``rubem -c `` spelling still works for one minor release and emits a ``DeprecationWarning``. @@ -21,7 +22,7 @@ from pydantic import ValidationError from . import __release__ -from ._deps import require_runtime_deps +from ._deps import require_calibration_deps, require_runtime_deps from .configuration._problems import ConfigurationError from .configuration.app_settings import AppSettings from .preprocessing.cli import app as preprocess_app @@ -178,6 +179,113 @@ def run( logger.info("RUBEM successfully finished!") +@app.command() +def calibrate( + configfile: Annotated[ + Path, + typer.Option( + "-c", + "--configfile", + callback=_configfile_callback, + help="Path to the configuration file (JSON).", + ), + ], + observed: Annotated[ + Path, + typer.Option( + "--observed", + exists=True, + dir_okay=False, + help="Observed series at the sample stations (CSV or PCRaster TSS).", + ), + ], + run_dir: Annotated[ + Path, + typer.Option("-o", "--run-dir", help="Directory for the calibration artifacts."), + ], + variable: Annotated[ + str, typer.Option("--variable", help="Output variable compared with the observed series.") + ] = "arn", + spinup_steps: Annotated[ + int, + typer.Option("--spinup-steps", min=0, help="Leading time steps excluded from the NSE."), + ] = 0, + maxiter: Annotated[ + int, typer.Option("--maxiter", min=1, help="Maximum number of generations.") + ] = 100, + popsize: Annotated[ + int, typer.Option("--popsize", min=1, help="Population multiplier of the search.") + ] = 15, + seed: Annotated[ + int | None, typer.Option("--seed", help="Seed of the differential evolution.") + ] = None, + workers: Annotated[ + int | None, + typer.Option("--workers", min=1, help="Parallel evaluations (default: all cores but one)."), + ] = None, + temp_dir: Annotated[ + Path | None, + typer.Option("--temp-dir", help="Parent of the per-evaluation output directories."), + ] = None, +) -> None: + """Calibrate the model parameters against an observed series.""" + require_runtime_deps() + require_calibration_deps() + + from .calibration.runner import CalibrationError, CalibrationSettings + from .calibration.runner import calibrate as run_calibration + + settings_arguments: dict = { + "variable": variable, + "spinup_steps": spinup_steps, + "maxiter": maxiter, + "popsize": popsize, + "seed": seed, + } + if workers is not None: + settings_arguments["workers"] = workers + if temp_dir is not None: + settings_arguments["temp_dir"] = str(temp_dir) + + try: + # The library only logs; the progress a command-line calibration is + # expected to show is written here, so that an embedded calibration + # stays silent unless its host configures logging. + print("Loading configuration and validating inputs...", flush=True) + print("Calibration started...", flush=True) + result = run_calibration( + configfile, observed, run_dir, CalibrationSettings(**settings_arguments) + ) + except (ConfigurationError, ValidationError, ValueError) as e: + # A configuration the user can fix: no traceback, the message says what. + logger.critical("Invalid configuration: %s", e) + raise typer.Exit(code=1) from e + except CalibrationError as e: + # The calibration itself refused to run or produced nothing usable; + # the message says what, so no traceback either. + logger.critical("Calibration failed: %s", e) + raise typer.Exit(code=1) from e + except KeyboardInterrupt as e: + logger.critical("RUBEM was interrupted by the user.") + raise typer.Exit(code=2) from e + except Exception as e: + logger.critical("RUBEM unexpectedly quit.") + logger.exception(e) + raise typer.Exit(code=1) from e + + best_nse = "n/a" if result.best_nse is None else f"{result.best_nse:.6f}" + print("Calibration finished successfully!", flush=True) + print(f"Best NSE: {best_nse}", flush=True) + print(f"Best objective: {result.best_objective:.6g}", flush=True) + print(f"Evaluations: {result.evaluations}", flush=True) + print(f"Generations: {result.generations}", flush=True) + print(f"Evaluations table: {result.evaluations_csv}", flush=True) + print(f"Result: {result.result_json}", flush=True) + print(f"Calibrated configuration: {result.calibrated_config}", flush=True) + + logger.info("RUBEM successfully finished!") + + @config_app.command("schema") def config_schema( schema_format: Annotated[ @@ -258,7 +366,7 @@ def _legacy_arguments(argv: Sequence[str]) -> list[str]: global _legacy_invocation _legacy_invocation = False arguments = list(argv) - if not arguments or arguments[0] in ("run", "config"): + if not arguments or arguments[0] in ("run", "config", "calibrate", "preprocess"): return arguments options = [a for a in arguments if a.startswith("-")] if not options or not {o.split("=", 1)[0] for o in options} <= _LEGACY_OPTIONS: diff --git a/tests/unit/test_cli_calibrate.py b/tests/unit/test_cli_calibrate.py new file mode 100644 index 00000000..3fb9f4f0 --- /dev/null +++ b/tests/unit/test_cli_calibrate.py @@ -0,0 +1,479 @@ +import json +from pathlib import Path + +import pytest + +from rubem.cli import main +from tests.helpers.synthetic import write_synthetic_dataset + + +@pytest.fixture(name="config_path") +def config_path_fixture(tmp_path): + path = tmp_path / "config.json" + path.write_text(json.dumps(write_synthetic_dataset(str(tmp_path))), encoding="utf8") + return path + + +@pytest.fixture(name="observed_path") +def observed_path_fixture(tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;1\n1;2.0\n2;3.0\n", encoding="utf8") + return path + + +class FakeCalibration: + """Stand in for ``rubem.calibration.runner.calibrate``. + + It records the arguments the command line built and answers with a result + whose files do not have to exist: the command only prints their paths. + """ + + def __init__(self, run_dir, best_nse=0.87): + from rubem.calibration.runner import CalibrationResult + + self.calls = [] + self.result = CalibrationResult( + best_parameters={ + "alpha": 4.5, + "beta": 0.5, + "w_1": 0.333, + "w_2": 0.333, + "w_3": 0.334, + "rcd": 5.0, + "f": 0.5, + "alpha_gw": 0.5, + "x": 0.5, + }, + best_nse=best_nse, + best_objective=169000.0, + evaluations=128, + generations=2, + success=True, + message="Optimization terminated successfully.", + run_dir=Path(run_dir), + evaluations_csv=Path(run_dir) / "evaluations.csv", + result_json=Path(run_dir) / "result.json", + calibrated_config=Path(run_dir) / "config-calibrated.json", + ) + + def __call__(self, config_path, observed_path, run_dir, settings): + self.calls.append((config_path, observed_path, run_dir, settings)) + return self.result + + +@pytest.fixture(name="fake_calibration") +def fake_calibration_fixture(tmp_path, monkeypatch): + fake = FakeCalibration(tmp_path / "calibration") + monkeypatch.setattr("rubem.calibration.runner.calibrate", fake) + return fake + + +class TestCliCalibrateArgumentsReachTheRunner: + @pytest.mark.unit + def test_the_paths_and_the_documented_defaults( + self, config_path, observed_path, tmp_path, fake_calibration, restore_logging + ): + run_dir = tmp_path / "calibration" + + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(run_dir), + ] + ) + + (given_config, given_observed, given_run_dir, settings) = fake_calibration.calls[0] + assert Path(given_config) == config_path + assert Path(given_observed) == observed_path + assert Path(given_run_dir) == run_dir + assert settings.variable == "arn" + assert settings.spinup_steps == 0 + assert settings.maxiter == 100 + assert settings.popsize == 15 + assert settings.seed is None + assert settings.temp_dir is None + + @pytest.mark.unit + def test_every_option_reaches_the_settings( + self, config_path, observed_path, tmp_path, fake_calibration, restore_logging + ): + temp_dir = tmp_path / "scratch" + + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + "--variable", + "rnf", + "--spinup-steps", + "12", + "--maxiter", + "3", + "--popsize", + "4", + "--seed", + "7", + "--workers", + "2", + "--temp-dir", + str(temp_dir), + ] + ) + + settings = fake_calibration.calls[0][3] + assert settings.variable == "rnf" + assert settings.spinup_steps == 12 + assert settings.maxiter == 3 + assert settings.popsize == 4 + assert settings.seed == 7 + assert settings.workers == 2 + assert Path(settings.temp_dir) == temp_dir + + @pytest.mark.unit + def test_the_runner_keeps_its_own_default_worker_count( + self, config_path, observed_path, tmp_path, fake_calibration, monkeypatch, restore_logging + ): + """Without ``--workers`` the command must not override the runner's default. + + The settings class is replaced by one whose default differs from the + ``None`` the option carries, so forwarding the option unconditionally + would be visible here. + """ + from dataclasses import dataclass + + from rubem.calibration.runner import CalibrationSettings + + @dataclass(frozen=True) + class SettingsWithAWorkerDefault(CalibrationSettings): + workers: int | None = 7 + + monkeypatch.setattr( + "rubem.calibration.runner.CalibrationSettings", SettingsWithAWorkerDefault + ) + + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + assert fake_calibration.calls[0][3].workers == 7 + + +class TestCliCalibrateOutput: + @pytest.mark.unit + def test_the_documented_progress_and_summary_are_printed( + self, config_path, observed_path, tmp_path, fake_calibration, capsys, restore_logging + ): + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + output = capsys.readouterr().out + result = fake_calibration.result + for expected in ( + "Loading configuration and validating inputs...", + "Calibration started...", + "Calibration finished successfully!", + "Best NSE: 0.87", + "Best objective: 169000", + "Evaluations: 128", + "Generations: 2", + str(result.evaluations_csv), + str(result.result_json), + str(result.calibrated_config), + ): + assert expected in output, f"missing line: {expected!r}\n{output}" + + @pytest.mark.unit + def test_a_run_without_a_usable_nse_says_so( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + fake = FakeCalibration(tmp_path / "calibration", best_nse=None) + monkeypatch.setattr("rubem.calibration.runner.calibrate", fake) + + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + assert "Best NSE: n/a" in capsys.readouterr().out + + +class TestCliCalibrateFailures: + @pytest.mark.unit + def test_a_calibration_error_exits_with_one_without_a_traceback( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + from rubem.calibration.runner import CalibrationError + + def fail(*args, **kwargs): + raise CalibrationError("every evaluation failed: no such variable") + + monkeypatch.setattr("rubem.calibration.runner.calibrate", fail) + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + captured = capsys.readouterr().err + assert error.value.code == 1 + assert "every evaluation failed: no such variable" in captured + assert "Traceback" not in captured + + @pytest.mark.unit + def test_an_unexpected_failure_logs_one_traceback( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + def fail(*args, **kwargs): + raise RuntimeError("the pool died") + + monkeypatch.setattr("rubem.calibration.runner.calibrate", fail) + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + captured = capsys.readouterr().err + assert error.value.code == 1 + assert captured.count("Traceback (most recent call last)") == 1 + assert "RUBEM unexpectedly quit." in captured + + @pytest.mark.unit + def test_a_missing_scipy_exits_with_the_installation_guidance( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + monkeypatch.setattr("rubem._deps.missing_calibration_deps", lambda: ["scipy"]) + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + message = str(error.value) + assert "scipy" in message + assert 'pip install "rubem[calibration]"' in message + + @pytest.mark.unit + def test_an_interrupted_calibration_exits_with_two( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + def interrupt(*args, **kwargs): + raise KeyboardInterrupt + + monkeypatch.setattr("rubem.calibration.runner.calibrate", interrupt) + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + captured = capsys.readouterr().err + assert error.value.code == 2 + assert "RUBEM was interrupted by the user." in captured + assert "Traceback" not in captured + + @pytest.mark.unit + @pytest.mark.parametrize( + ("option", "value"), + [ + ("--maxiter", "0"), + ("--maxiter", "-1"), + ("--popsize", "0"), + ("--popsize", "-1"), + ("--workers", "0"), + ("--workers", "-1"), + ("--spinup-steps", "-1"), + ], + ) + def test_a_search_size_outside_its_range_is_rejected( + self, config_path, observed_path, tmp_path, fake_calibration, option, value, restore_logging + ): + """A generation count, a population and a worker count all start at one.""" + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + option, + value, + ] + ) + + assert error.value.code == 2 + assert not fake_calibration.calls + + @pytest.mark.unit + def test_a_negative_spin_up_is_rejected( + self, config_path, observed_path, tmp_path, fake_calibration, capsys, restore_logging + ): + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + "--spinup-steps", + "-1", + ] + ) + + assert error.value.code == 2 + assert not fake_calibration.calls + + @pytest.mark.unit + def test_a_missing_observed_file_is_rejected_before_anything_runs( + self, config_path, tmp_path, fake_calibration, capsys, restore_logging + ): + missing = tmp_path / "missing.csv" + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(missing), + "-o", + str(tmp_path / "calibration"), + ] + ) + + assert error.value.code == 2 + assert not fake_calibration.calls + assert "Loading configuration" not in capsys.readouterr().out + + @pytest.mark.unit + def test_an_incomplete_command_is_not_mapped_onto_the_legacy_spelling( + self, config_path, capsys, restore_logging + ): + """``calibrate`` is a subcommand, not the deprecated bare ``-c`` form.""" + with pytest.raises(SystemExit) as error: + main(["calibrate", "-c", str(config_path)]) + + captured = capsys.readouterr() + assert error.value.code == 2 + assert "Missing option '--observed'" in captured.err + assert "is deprecated" not in captured.err + + +class TestCliCalibrateHelp: + @pytest.mark.unit + def test_the_command_is_listed_in_the_help(self, capsys, restore_logging): + with pytest.raises(SystemExit) as error: + main(["--help"]) + + assert error.value.code == 0 + assert "calibrate" in capsys.readouterr().out + + @pytest.mark.unit + def test_the_help_documents_every_option(self, capsys, restore_logging): + with pytest.raises(SystemExit) as error: + main(["calibrate", "--help"]) + + assert error.value.code == 0 + output = capsys.readouterr().out + assert "Calibrate the model parameters against an observed series." in output + for option in ( + "--configfile", + "--observed", + "--run-dir", + "--variable", + "--spinup-steps", + "--maxiter", + "--popsize", + "--seed", + "--workers", + "--temp-dir", + ): + assert option in output, f"missing option: {option}\n{output}" + + @pytest.mark.unit + def test_importing_the_command_line_does_not_need_scipy(self): + """The CLI module stays importable without the optional extra.""" + import subprocess + import sys + + completed = subprocess.run( + [sys.executable, "-c", "import sys, rubem.cli; print('scipy' in sys.modules)"], + cwd=Path(__file__).resolve().parents[2], + capture_output=True, + text=True, + check=True, + ) + assert completed.stdout.strip() == "False" diff --git a/tests/unit/test_deps.py b/tests/unit/test_deps.py index 410d5bca..41bcb5ae 100644 --- a/tests/unit/test_deps.py +++ b/tests/unit/test_deps.py @@ -18,3 +18,54 @@ def test_missing_dependencies_are_named_with_an_actionable_source(self, monkeypa assert name in message assert "conda-forge" in message assert _deps._ENVIRONMENT_YML_URL in message + + +class TestRequireCalibrationDeps: + @pytest.mark.unit + def test_present_dependencies_do_not_raise(self): + _deps.require_calibration_deps() + + @pytest.mark.unit + def test_nothing_is_missing_when_every_dependency_is_importable(self): + assert _deps.missing_calibration_deps() == [] + + @pytest.mark.unit + def test_an_unimportable_dependency_is_reported_as_missing(self, monkeypatch): + monkeypatch.setattr(_deps.importlib.util, "find_spec", lambda name: None) + + assert _deps.missing_calibration_deps() == list(_deps._CALIBRATION_DEPENDENCIES) + + @pytest.mark.unit + def test_a_broken_dependency_counts_as_missing(self, monkeypatch): + def raise_import_error(name): + raise ImportError(name) + + monkeypatch.setattr(_deps.importlib.util, "find_spec", raise_import_error) + + assert _deps.missing_calibration_deps() == list(_deps._CALIBRATION_DEPENDENCIES) + + @pytest.mark.unit + def test_the_message_names_the_dependency_and_the_extra(self): + message = _deps.calibration_deps_message(["scipy"]) + + assert "scipy" in message + assert 'pip install "rubem[calibration]"' in message + + @pytest.mark.unit + def test_the_message_defaults_to_every_calibration_dependency(self): + """The runner asks for the guidance without having computed the list.""" + assert _deps.calibration_deps_message() == _deps.calibration_deps_message( + list(_deps._CALIBRATION_DEPENDENCIES) + ) + + @pytest.mark.unit + def test_missing_dependencies_are_named_with_an_actionable_extra(self, monkeypatch): + monkeypatch.setattr(_deps.importlib.util, "find_spec", lambda name: None) + + with pytest.raises(SystemExit) as error: + _deps.require_calibration_deps() + + message = str(error.value) + for name in _deps._CALIBRATION_DEPENDENCIES: + assert name in message + assert 'pip install "rubem[calibration]"' in message From cdd4a4188c819471c386003fef4f1fbc15748720 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Wed, 16 Sep 2026 17:50:40 -0300 Subject: [PATCH 04/12] Package SciPy as the optional calibration extra Add the rubem[calibration] extra (scipy>=1.16.2, the first release with CPython 3.14 wheels on Linux and Windows), include it in the dev extra and in the CI test installation, commit the ci/dev-requirements.in the audit lock was compiled from, with scipy added, and recompile the lock keeping every existing pin. --- .github/workflows/ci.yml | 2 +- ci/dev-requirements.in | 18 ++++++++++++++++++ ci/dev-requirements.lock | 34 ++++++++++++++++++---------------- pyproject.toml | 7 ++++++- 4 files changed, 43 insertions(+), 18 deletions(-) create mode 100644 ci/dev-requirements.in diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1366150b..53b27abd 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -113,7 +113,7 @@ jobs: MATRIX_PYTHON: ${{ matrix.python }} shell: bash -el {0} - name: Install RUBEM with test dependencies - run: pip install -e ".[test,docs]" + run: pip install -e ".[test,docs,calibration]" shell: bash -el {0} - name: Run the test suite run: pytest -n auto -ra --cov=rubem --cov-report=xml --junitxml="TestResults-Rubem.xml" diff --git a/ci/dev-requirements.in b/ci/dev-requirements.in new file mode 100644 index 00000000..eec9fe2e --- /dev/null +++ b/ci/dev-requirements.in @@ -0,0 +1,18 @@ +# Direct development requirements audited by the pip-audit job. Refresh the +# lock with: +# uv pip compile ci/dev-requirements.in -o ci/dev-requirements.lock --python-version 3.12 +# An existing lock keeps its pins; add --upgrade to move every package forward. +build +codespell +pre-commit +pyfakefs +pytest-cov +pytest-mock +pytest-xdist +ruff +scipy>=1.16.2 +sphinx-autodoc-typehints +sphinx-rtd-theme +sphinxcontrib-email +sphinxcontrib-youtube +twine diff --git a/ci/dev-requirements.lock b/ci/dev-requirements.lock index 0e3042c7..00e304db 100644 --- a/ci/dev-requirements.lock +++ b/ci/dev-requirements.lock @@ -1,11 +1,11 @@ # This file was autogenerated by uv via the following command: -# uv pip compile dev-requirements.in -o dev-requirements.lock --python-version 3.12 +# uv pip compile ci/dev-requirements.in -o ci/dev-requirements.lock --python-version 3.12 alabaster==1.0.0 # via sphinx babel==2.18.0 # via sphinx build==1.5.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in certifi==2026.7.22 # via requests cffi==2.1.1 @@ -15,7 +15,7 @@ cfgv==3.5.0 charset-normalizer==3.5.1 # via requests codespell==2.4.3 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in coverage==7.15.4 # via pytest-cov cryptography==50.0.0 @@ -73,6 +73,8 @@ nh3==0.3.7 # via readme-renderer nodeenv==1.10.0 # via pre-commit +numpy==2.5.3 + # via scipy packaging==26.3 # via # build @@ -86,11 +88,11 @@ pluggy==1.6.0 # pytest # pytest-cov pre-commit==4.6.2 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in pycparser==3.0 # via cffi pyfakefs==6.2.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in pygments==2.21.0 # via # pytest @@ -101,16 +103,15 @@ pyproject-hooks==1.2.0 # via build pytest==9.1.1 # via - # -r dev-requirements.in # pytest-cov # pytest-mock # pytest-xdist pytest-cov==7.1.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in pytest-mock==3.15.1 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in pytest-xdist==3.8.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in python-discovery==1.5.2 # via virtualenv pyyaml==6.0.3 @@ -132,29 +133,30 @@ rich==15.0.0 roman-numerals==4.1.0 # via sphinx ruff==0.16.4 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in +scipy==1.18.1 + # via -r ci/dev-requirements.in secretstorage==3.5.0 # via keyring snowballstemmer==3.1.1 # via sphinx sphinx==9.1.0 # via - # -r dev-requirements.in # sphinx-autodoc-typehints # sphinx-rtd-theme # sphinxcontrib-email # sphinxcontrib-jquery # sphinxcontrib-youtube sphinx-autodoc-typehints==3.13.2 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in sphinx-rtd-theme==3.1.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in sphinxcontrib-applehelp==2.0.0 # via sphinx sphinxcontrib-devhelp==2.0.0 # via sphinx sphinxcontrib-email==0.3.6 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in sphinxcontrib-htmlhelp==2.1.0 # via sphinx sphinxcontrib-jquery==4.1 @@ -166,9 +168,9 @@ sphinxcontrib-qthelp==2.0.0 sphinxcontrib-serializinghtml==2.0.0 # via sphinx sphinxcontrib-youtube==1.5.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in twine==7.0.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in urllib3==2.7.0 # via # id diff --git a/pyproject.toml b/pyproject.toml index 14df2be4..be714733 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -66,8 +66,13 @@ preprocessing = [ "pykrige>=1.7", "scikit-gstat>=1.0", ] +# 1.16.2 is the first release with CPython 3.14 wheels for Linux and Windows, +# the range of interpreters this project supports. +calibration = [ + "scipy>=1.16.2", +] dev = [ - "rubem[test,docs,preprocessing]", + "rubem[test,docs,preprocessing,calibration]", "build", "codespell", "pre-commit", From b7ace36151cfc6c330e6d7ff824a4fe421ca16e8 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Wed, 16 Sep 2026 17:50:40 -0300 Subject: [PATCH 05/12] Document the calibrate command in the user guide --- doc/source/userguide.rst | 62 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 57 insertions(+), 5 deletions(-) diff --git a/doc/source/userguide.rst b/doc/source/userguide.rst index aa39b9d1..b53ea134 100644 --- a/doc/source/userguide.rst +++ b/doc/source/userguide.rst @@ -1127,9 +1127,11 @@ Configuration File Template Running RUBEM ------------- -The ``rubem`` command has two subcommands: ``run`` executes a simulation and -``config schema`` prints the JSON Schema of the configuration file. Running it -without arguments shows the available commands: +The ``rubem`` command groups its work in subcommands: ``run`` executes a +simulation, ``calibrate`` fits the calibration parameters to an observed +series, ``config schema`` prints the JSON Schema of the configuration file and +``preprocess`` prepares the input rasters. Running it without arguments shows +the available commands: .. code-block:: console @@ -1143,8 +1145,10 @@ without arguments shows the available commands: -h, --help Show this message and exit. Commands: - run Run a simulation from a configuration file. - config Inspect the configuration file format. + run Run a simulation from a configuration file. + calibrate Calibrate the model parameters against an observed series. + config Inspect the configuration file format. + preprocess Prepare model inputs: inspect and convert rasters, build... .. note:: @@ -1170,6 +1174,54 @@ Use ``-h`` or ``--help`` to get a brief description of each command and its argu the model. -h, --help Show this message and exit. +.. code-block:: console + + $ rubem calibrate -h + Usage: rubem calibrate [OPTIONS] + + Calibrate the model parameters against an observed series. + + Options: + -c, --configfile Path to the configuration file (JSON). + [required] + --observed Observed series at the sample stations (CSV or + PCRaster TSS). [required] + -o, --run-dir Directory for the calibration artifacts. + [required] + --variable Output variable compared with the observed + series. [default: arn] + --spinup-steps Leading time steps excluded from the NSE. + [default: 0; x>=0] + --maxiter Maximum number of generations. [default: 100; + x>=1] + --popsize Population multiplier of the search. [default: + 15; x>=1] + --seed Seed of the differential evolution. + --workers Parallel evaluations (default: all cores but + one). [x>=1] + --temp-dir Parent of the per-evaluation output directories. + -h, --help Show this message and exit. + +``calibrate`` needs SciPy, which comes with ``pip install +"rubem[calibration]"``. It writes ``evaluations.csv``, ``result.json`` +and ``-calibrated.json`` into the ``--run-dir``. The method itself (the +objective, the parameter bounds, the derivation of ``w3`` and how to budget the +number of evaluations) is documented separately. + +.. code-block:: console + + $ rubem calibrate -c project-config.json --observed observed.csv -o calibration + Loading configuration and validating inputs... + Calibration started... + Calibration finished successfully! + Best NSE: 0.873210 + Best objective: 160757 + Evaluations: 2048 + Generations: 16 + Evaluations table: /basins/ipojuca/calibration/evaluations.csv + Result: /basins/ipojuca/calibration/result.json + Calibrated configuration: /basins/ipojuca/calibration/project-config-calibrated.json + Use ``-V`` or ``--version`` to get the version of the RUBEM. .. code-block:: console From a99e1049ee26e9e818f04b266a676d35a57b5c23 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Wed, 16 Sep 2026 17:50:40 -0300 Subject: [PATCH 06/12] Record the calibrate command in the changelog --- doc/source/changelog.rst | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index e39985a1..385dccf9 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -13,6 +13,15 @@ The format follows `Keep a Changelog `__. Added ````` +- Added ``rubem calibrate``: a differential evolution (SciPy) over the eight free + calibration parameters, ``w3`` derived from the other two weights, minimizing + ``1000 (100 (1 - NSE))^2`` on the station series of one output variable + (``arn`` by default) against an observed CSV or PCRaster TSS file, with + masked values, a zero-variance guard and a spin-up window; every candidate + runs in a fresh spawned worker through the Python API, records one JSON + file, and the run directory receives ``evaluations.csv``, ``result.json`` + and ``-calibrated.json``. SciPy comes with the optional extra + ``rubem[calibration]`` (`#343 `__). - Added ``rubem.api``, the public Python surface: ``Model.from_file`` and ``Model.from_config`` load a configuration file or document, ``Model.run()`` runs the simulation in the current process and ``Model.run_isolated()`` in a From deb85e2bcaa6803cb5772d3113352b56f6c84ca2 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Wed, 16 Sep 2026 18:49:54 -0300 Subject: [PATCH 07/12] Refuse an observed series that matches no station or step before the search An observed series that shares no time step with the simulated steps after the spin-up, or no station id with the sample locations raster, makes every evaluation fail after its simulation has run, so a wrong file cost a whole generation of model runs before the calibration failed. The parent now settles both cases before the run directory is prepared and before any worker starts, warns once about observed stations the configuration does not sample, and skips the station check for the zones aggregation, whose ids are only known at run time. The reproducibility statement of the runner is qualified: ordering the rows by candidate removes the scheduling order, but two runs of one seed evaluate the same candidates only when the simulation itself is reproducible, which on the real basins needs a fixed LDD raster. --- rubem/calibration/runner.py | 155 ++++++++++++++++++++++-- tests/unit/calibration/test_runner.py | 165 ++++++++++++++++++++++++-- 2 files changed, 304 insertions(+), 16 deletions(-) diff --git a/rubem/calibration/runner.py b/rubem/calibration/runner.py index 4f99fb6f..3d73fd21 100644 --- a/rubem/calibration/runner.py +++ b/rubem/calibration/runner.py @@ -11,8 +11,15 @@ evaluation, ``result.json``, the best candidate and the settings that found it, and ``-calibrated.json``, the calibrated configuration in the format the input was written in. The rows of the table are ordered by the -candidate they evaluated and not by the moment they were written, so that two -calibrations of one configuration with one seed produce the same table. +candidate they evaluated and not by the moment they were written, so that the +order of the rows does not depend on the order in which the workers happened to +finish. It does not make the table itself reproducible: the ``id``, ``pid`` and +``elapsed_seconds`` columns differ between two runs of one seed, and from the +first generation on the selection reads the objective values, so two runs +evaluate the same candidates only when the simulation itself is reproducible +(on the real basins that requires a fixed ``RASTERS.ldd``, since ``lddcreate`` +does not always derive the same directions twice; see the note on the LDD +raster in the user guide). Number of evaluations The population of a ``sobol`` initialization is ``max(5, popsize * 8)`` @@ -51,7 +58,7 @@ ModelConfigurationFileV1, ) from ._worker import _INADMISSIBLE_ERROR, EvaluationContext, evaluate -from .objective import read_series +from .objective import Series, read_series from .parameters import ( CALIBRATION_PARAMETERS, FREE_PARAMETERS, @@ -293,7 +300,9 @@ def calibrate( :rtype: CalibrationResult :raises CalibrationError: If SciPy is not installed, if the configuration - cannot be calibrated as asked, or if every evaluation of the run failed. + cannot be calibrated as asked, if the observed series shares no time + step or no station with what the configuration will sample, or if every + evaluation of the run failed. :raises ImportError: If PCRaster or GDAL are not installed. :raises FileNotFoundError: If the configuration or the observed series is not there. @@ -323,6 +332,19 @@ def calibrate( _check_sample_locations(configuration, file_v1) observed = read_series(observed_path) + aggregation = file_v1.model_simulation_output.time_series_samples.aggregation + _check_observed_series( + observed, + first_step=configuration.simulation_period.first_step, + last_step=configuration.simulation_period.last_step, + spinup_steps=settings.spinup_steps, + aggregation=aggregation, + station_ids=( + None + if aggregation is Aggregation.ZONES + else _station_ids(configuration.raster_files.sample_locations) + ), + ) x0 = parameters_to_vector(configuration.calibration_parameters.model_dump()) run_directory = as_path(run_dir).absolute() @@ -500,6 +522,118 @@ def _check_sample_locations(configuration, file_v1: ModelConfigurationFileV1) -> ) +def _station_ids(sample_locations: PathInput) -> set[str]: + """Return the station ids of a sample locations raster. + + The raster is a nominal map, so its cells read as integers. Neither the + no-data cells nor the zeroes are stations, which is what the validation of + the inputs already asks of a sample raster: its ids are ``1..N`` and its + background is missing or zero. Every remaining value becomes the string the + model writes in the header of its time series table, ``str(int(value))``. + + :param sample_locations: The sample locations raster of the configuration. + :type sample_locations: str | os.PathLike[str] + + :return: The distinct station ids, as strings. + :rtype: set[str] + """ + from osgeo import gdal + + gdal.UseExceptions() + gdal.AllRegister() + dataset = gdal.Open(str(as_path(sample_locations))) + try: + band = dataset.GetRasterBand(1) + values = np.unique(np.asarray(band.ReadAsArray())) + no_data = band.GetNoDataValue() + finally: + dataset = None + + if no_data is not None: + values = values[values != no_data] + return {str(int(value)) for value in values.tolist() if math.isfinite(value) and value != 0} + + +def _check_observed_series( + observed: Series, + *, + first_step: int, + last_step: int, + spinup_steps: int, + aggregation: Aggregation, + station_ids: set[str] | None, +) -> None: + """Refuse an observed series that cannot be compared with what the run samples. + + Every evaluation aligns the observed series with the table its own + simulation wrote, so a series that shares no step or no station with that + table makes every evaluation fail. The mismatch is the same for every + candidate, so it is settled here, in the parent, before a single model run + is spent on it. + + The compared steps are the simulated steps, ``first_step`` to ``last_step``, + that survive the spin-up. The stations are the ids the ``point`` and the + ``subcatchment`` aggregations sample, the values of the sample locations + raster; the ``zones`` aggregation remaps the zone ids to ``1..N`` at run + time (the mapping is written to ``zones_mapping.csv``), so the parent cannot + know the ids of its stations and the station check is skipped for it. + + :param observed: The observed series, as the parent read it. + :type observed: rubem.calibration.objective.Series + + :param first_step: Number of the first simulated step. + :type first_step: int + + :param last_step: Number of the last simulated step. + :type last_step: int + + :param spinup_steps: Number of initial time steps excluded from the + efficiency. + :type spinup_steps: int + + :param aggregation: The aggregation the time series of the run use. + :type aggregation: Aggregation + + :param station_ids: The ids of the stations the run will sample, or ``None`` + when they are not known before the run. + :type station_ids: set[str] | None + + :raises CalibrationError: If no observed step is compared, or if no observed + station id is one of the stations of the run. + """ + compared_steps = {step for step in range(first_step, last_step + 1) if step > spinup_steps} + observed_steps = {int(step) for step in np.asarray(observed.steps).tolist()} + if not compared_steps & observed_steps: + raise CalibrationError( + f"The observed series covers the step(s) {min(observed_steps)} to " + f"{max(observed_steps)}, and the configuration simulates the step(s) " + f"{first_step} to {last_step}, of which the {spinup_steps} spin-up step(s) " + "are excluded from the efficiency: the two series would have no time step " + "in common and every evaluation would fail." + ) + + if aggregation is Aggregation.ZONES or station_ids is None: + return + + observed_ids = set(observed.stations) + if not observed_ids & station_ids: + raise CalibrationError( + f"The observed series has the station(s) {', '.join(sorted(observed_ids)) or 'none'} " + f"and the configuration samples the station(s) " + f"{', '.join(sorted(station_ids)) or 'none'}: the two series would have no " + "station in common and every evaluation would fail." + ) + + foreign = sorted(observed_ids - station_ids) + if foreign: + logger.warning( + "The observed series has %d station(s) the configuration does not sample " + "(%s); they are ignored by the calibration.", + len(foreign), + ", ".join(foreign), + ) + + def _population_size(settings: CalibrationSettings) -> int: """Return the number of population members the search will use. @@ -541,10 +675,15 @@ def _record_order(record: dict[str, Any]) -> tuple[tuple[float, ...], str]: The records are not ordered by the moment they were written: the workers run in parallel, so that order changes from run to run. Ordering them by the - candidate they evaluated makes ``evaluations.csv`` a function of the seed - alone, so two calibrations of the same configuration with the same seed - produce the same table and can be compared line by line. The identifier only - separates two evaluations of the very same candidate. + candidate they evaluated takes that scheduling order out of + ``evaluations.csv``. It does not make the table reproducible on its own: the + ``id``, ``pid`` and ``elapsed_seconds`` columns differ between two runs, and + from the first generation on the selection reads the objective values, so + two runs of one seed evaluate the same candidates only when the simulation + itself is reproducible (on the real basins that requires a fixed + ``RASTERS.ldd``, since ``lddcreate`` does not always derive the same + directions twice; see the note on the LDD raster in the user guide). The + identifier only separates two evaluations of the very same candidate. """ parameters = record.get("parameters") or {} values = tuple(float(parameters.get(name, math.inf)) for name in CALIBRATION_PARAMETERS) diff --git a/tests/unit/calibration/test_runner.py b/tests/unit/calibration/test_runner.py index d1fe1b53..0d340619 100644 --- a/tests/unit/calibration/test_runner.py +++ b/tests/unit/calibration/test_runner.py @@ -18,6 +18,7 @@ from rubem import _deps from rubem.api import Model from rubem.calibration import runner +from rubem.calibration.objective import Series from rubem.calibration.parameters import ( CALIBRATION_PARAMETERS, FREE_PARAMETERS, @@ -28,13 +29,18 @@ EVALUATION_COLUMNS, CalibrationError, CalibrationSettings, + _check_observed_series, _population_size, _Progress, + _station_ids, calibrate, ) from rubem.configuration.model_configuration import ModelConfiguration from rubem.configuration.model_configuration_file import ModelConfigurationFile -from rubem.configuration.model_configuration_file_v1 import ModelConfigurationFileV1 +from rubem.configuration.model_configuration_file_v1 import ( + Aggregation, + ModelConfigurationFileV1, +) from tests.helpers.config import REPO_ROOT from tests.helpers.synthetic import write_synthetic_dataset @@ -362,9 +368,10 @@ def candidates(result): for row in read_evaluations(result.evaluations_csv) ] - # Same seed, same initial population, same candidates, in the same order: - # the table of a calibration is a function of its settings alone. Only - # the identifiers and the process ids differ between the two runs. + # Same seed, same initial population, same candidates, in the same + # order: the synthetic dataset carries a fixed LDD, so the simulation + # itself is reproducible and the two searches follow the same path. Only + # the identifiers, the process ids and the elapsed times differ. assert candidates(first) == candidates(second) assert first.best_parameters == second.best_parameters assert first.best_objective == second.best_objective @@ -385,11 +392,21 @@ class TestFailures: def test_a_run_whose_evaluations_all_fail_is_reported_with_its_first_error( self, dataset, tmp_path ): - foreign = tmp_path / "foreign.csv" - foreign.write_text("0;A;B\n1;1.0;2.0\n2;2.0;3.0\n3;3.0;4.0\n", encoding="utf8") + # The stations and the steps of this file are the ones the run samples, + # so the parent lets it through; every value is a gap, so every worker + # fails inside the evaluation of the series. + gaps = tmp_path / "gaps.csv" + gaps.write_text( + "0;1;2\n" + "".join(f"{step};-9999;-9999\n" for step in range(1, TIMESTEPS + 1)), + encoding="utf8", + ) - with pytest.raises(CalibrationError, match="no station in common"): - dataset.calibrate(observed=foreign, maxiter=1) + with pytest.raises(CalibrationError, match="Nash-Sutcliffe") as failure: + dataset.calibrate(observed=gaps, maxiter=1) + + # The search ran: what is reported is the summary of the finished run, + # carrying the first error a worker recorded, not a refusal of the parent. + assert "Every evaluation of the calibration failed" in str(failure.value) rows = read_evaluations(dataset.run_dir / "evaluations.csv") assert rows, "the table of the evaluations is written even when every one failed" @@ -430,6 +447,138 @@ def test_without_scipy_the_optional_extra_is_named(self, dataset, monkeypatch): assert not dataset.run_dir.exists(), "nothing is written when the search cannot start" +@pytest.fixture +def no_pool(monkeypatch): + """Let a test prove that no worker process was ever started.""" + _RecordingExecutor.built.clear() + monkeypatch.setattr(runner, "ProcessPoolExecutor", _RecordingExecutor) + return _RecordingExecutor.built + + +class TestObservedSeries: + @pytest.mark.unit + def test_an_observed_series_of_foreign_stations_is_refused_before_the_search( + self, dataset, tmp_path, no_pool + ): + foreign = tmp_path / "foreign.csv" + foreign.write_text("0;A;B\n1;1.0;2.0\n2;2.0;3.0\n3;3.0;4.0\n", encoding="utf8") + + with pytest.raises(CalibrationError, match="no station in common") as failure: + dataset.calibrate(observed=foreign, maxiter=1) + + # Both sets are named, and the refusal comes before anything is written + # or any worker is started: a wrong observed file costs no model run. + assert "A, B" in str(failure.value) + assert "1, 2" in str(failure.value) + assert not dataset.run_dir.exists() + assert no_pool == [] + assert multiprocessing.active_children() == [] + + @pytest.mark.unit + @pytest.mark.parametrize( + ("steps", "overrides"), + [ + ([TIMESTEPS + 1, TIMESTEPS + 2], {}), + ([-1, 0], {}), + ([1, TIMESTEPS], {"spinup_steps": TIMESTEPS}), + ], + ids=["after the window", "before the window", "inside the spin-up"], + ) + def test_an_observed_series_outside_the_compared_steps_is_refused( + self, dataset, tmp_path, no_pool, steps, overrides + ): + elsewhere = tmp_path / "elsewhere.csv" + elsewhere.write_text( + "0;1;2\n" + "".join(f"{step};1.0;2.0\n" for step in steps), encoding="utf8" + ) + + with pytest.raises(CalibrationError, match="no time step in common") as failure: + dataset.calibrate(observed=elsewhere, **overrides) + + message = str(failure.value) + assert f"{min(steps)} to {max(steps)}" in message + assert f"1 to {TIMESTEPS}" in message + assert f"{overrides.get('spinup_steps', 0)} spin-up" in message + assert not dataset.run_dir.exists() + assert no_pool == [] + + @pytest.mark.unit + def test_an_observed_station_the_run_does_not_sample_is_warned_about_and_ignored( + self, dataset, caplog + ): + rows = dataset.observed.read_text(encoding="utf8").splitlines() + extra = dataset.observed.with_name("extra.csv") + extra.write_text( + "".join(f"{row};{'3' if number == 0 else '1.0'}\n" for number, row in enumerate(rows)), + encoding="utf8", + ) + + exact = calibrate( + dataset.config_file, dataset.observed, dataset.run_dir / "exact", dataset.settings() + ) + with caplog.at_level(logging.WARNING, logger="rubem.calibration.runner"): + with_extra = calibrate( + dataset.config_file, extra, dataset.run_dir / "extra", dataset.settings() + ) + + # The foreign column is named once and then ignored by the alignment, so + # the calibration is the one the exact file produces. + assert "does not sample (3)" in caplog.text + assert with_extra.best_parameters == exact.best_parameters + assert with_extra.best_objective == exact.best_objective + assert with_extra.best_nse == exact.best_nse + assert with_extra.evaluations == exact.evaluations + + @pytest.mark.unit + def test_the_stations_are_the_values_of_the_sample_locations_raster(self, dataset): + samples = ModelConfiguration( + dataset.config, validate_input=False + ).raster_files.sample_locations + + # The no-data cells of the nominal map are not stations, and the ids are + # the strings the model writes in the header of its table. + assert _station_ids(samples) == {"1", "2"} + + @pytest.mark.unit + def test_a_background_of_zeroes_is_not_a_station(self, tmp_path): + from osgeo import gdal + + gdal.UseExceptions() + gdal.AllRegister() + raster = tmp_path / "samples-with-zeroes.tif" + handle = gdal.GetDriverByName("GTiff").Create(str(raster), 3, 3, 1, gdal.GDT_Int32) + try: + values = np.zeros((3, 3), dtype=np.int32) + values[0, 0], values[2, 2] = 1, 2 + handle.GetRasterBand(1).WriteArray(values) + handle.FlushCache() + finally: + handle = None + + # A raster may mark its background with zeroes and declare no no-data + # value at all; the validation of the inputs counts neither the zeroes + # nor the missing cells as a sample, and neither does this. + assert _station_ids(raster) == {"1", "2"} + + @pytest.mark.unit + def test_the_stations_of_a_zones_aggregation_are_not_checked(self): + # The model remaps the zone ids to 1..N at run time, so the parent has + # no station id to compare and the check only covers the steps. + observed = Series(steps=np.array([1, 2, 3]), stations={"A": np.zeros(3)}) + + assert ( + _check_observed_series( + observed, + first_step=1, + last_step=3, + spinup_steps=0, + aggregation=Aggregation.ZONES, + station_ids={"1", "2"}, + ) + is None + ) + + class TestSettings: @pytest.mark.unit def test_the_default_decision_vector_is_the_configuration_under_calibration(self, dataset): From 4f9387f551f486a93c010151bb629a45a782e6b1 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Wed, 16 Sep 2026 18:30:13 -0300 Subject: [PATCH 08/12] Add the dataset marker and the RUBEM_DATASET_DIR convention Tests marked dataset need the local copy of the published datasets, which the environment variable names; when it is unset, or does not name a directory, they are skipped with the reason, whether or not they request the dataset_dir fixture, so the suite stays self-contained elsewhere. --- pyproject.toml | 1 + tests/conftest.py | 60 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index be714733..29b24f6a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -95,6 +95,7 @@ markers = [ "unit: mark a test as a unit test.", "integration: mark a test as a integration test.", "exact: mark a test that reproduces the goldens byte-for-byte on the golden environment.", + "dataset: mark a test that needs the published datasets (RUBEM_DATASET_DIR).", ] [tool.coverage.run] diff --git a/tests/conftest.py b/tests/conftest.py index eb861b6f..cd305f0f 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,8 +1,68 @@ import logging import os +from pathlib import Path import pytest +DATASET_DIR_VARIABLE = "RUBEM_DATASET_DIR" +"""Environment variable pointing at the local copy of the published datasets. + +The directory it names holds the manifest ``_basin.json`` of every +published basin beside the extracted ``_basin/`` directory of its inputs. +The datasets are hundreds of megabytes and are not part of the repository, so +the tests that need them are marked ``dataset`` and run only where the variable +is set; everywhere else the suite stays self-contained. +""" + +NO_DATASET_DIR = f"{DATASET_DIR_VARIABLE} is not set" +"""The reason a test marked ``dataset`` is skipped with when the variable is unset.""" + + +def _dataset_dir() -> tuple[Path | None, str]: + """Return the datasets directory, or ``None`` and the reason it is not available. + + An unset or empty variable means the published datasets are not on this + machine; a variable that names something other than a directory is reported + as such, so that a mistyped path or a mount that is not up does not read as + "not requested". + """ + value = os.environ.get(DATASET_DIR_VARIABLE) + if not value: + return None, NO_DATASET_DIR + directory = Path(value) + if not directory.is_dir(): + return None, f"{DATASET_DIR_VARIABLE}={value} is not a directory" + return directory, "" + + +def pytest_collection_modifyitems(config, items): + """Skip every test marked ``dataset`` when the datasets are not there. + + The marker alone is enough: a test that never asks for the ``dataset_dir`` + fixture, because it reaches the datasets through another fixture or through + a helper, is skipped for the same reason as one that does. + """ + directory, reason = _dataset_dir() + if directory is not None: + return + skip = pytest.mark.skip(reason=reason) + for item in items: + if item.get_closest_marker("dataset") is not None: + item.add_marker(skip) + + +@pytest.fixture +def dataset_dir(): + """The directory of the published datasets, or a skipped test. + + :return: The directory :data:`DATASET_DIR_VARIABLE` names. + :rtype: pathlib.Path + """ + directory, reason = _dataset_dir() + if directory is None: + pytest.skip(reason) + return directory + @pytest.fixture(autouse=True) def preserve_working_directory(): From 9e2cbddfe2d7f00e80c5c7144a5e2eaa81a2fb60 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Wed, 16 Sep 2026 18:30:13 -0300 Subject: [PATCH 09/12] Add the calibration smoke test on the Ipojuca basin Rebase the manifest onto the local copy, restrict the window to six monthly steps, generate the drainage network once and name it in the configuration (the accumulated runoff is routed on a network lddcreate does not reproduce between calls), disable GDAL's auxiliary metadata so nothing is written beside the dataset, take the series of a base run as the observations, and check that a one-generation differential evolution from five members finds the configuration again with an objective of exactly zero, consistent artifacts and no process left behind. --- tests/dataset/__init__.py | 0 tests/dataset/test_calibrate_ipojuca.py | 200 ++++++++++++++++++++++++ 2 files changed, 200 insertions(+) create mode 100644 tests/dataset/__init__.py create mode 100644 tests/dataset/test_calibrate_ipojuca.py diff --git a/tests/dataset/__init__.py b/tests/dataset/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/dataset/test_calibrate_ipojuca.py b/tests/dataset/test_calibrate_ipojuca.py new file mode 100644 index 00000000..4f39c89c --- /dev/null +++ b/tests/dataset/test_calibrate_ipojuca.py @@ -0,0 +1,200 @@ +"""Smoke test of the calibration on the published Ipojuca basin dataset. + +The test rebuilds the configuration of the basin from its manifest, runs the +model once over a short window to obtain the series a perfect calibration has +to reproduce, and then asks the differential evolution to find the parameters +that produced it again, from a five-member initial population, in one +generation. + +The drainage network is fixed before anything runs. ``lddcreate`` is called by +the model whenever the configuration names no ``ldd`` raster, and on the real +basin digital elevation models it returns a different network on every call, so +the accumulated runoff, which is routed along it, is not reproducible between +two runs of one configuration. A calibration on ``arn`` without a fixed network +would compare series routed on different networks; the test therefore generates +the network once and names it in the configuration, as the documentation of the +command requires. + +The dataset is not part of the repository: the whole module is marked +``dataset`` and runs only where ``RUBEM_DATASET_DIR`` points at the local copy. +""" + +import csv +import json +import multiprocessing +import shutil +from pathlib import Path + +import numpy as np +import pytest + +from rubem.api import Model +from rubem.calibration.parameters import CALIBRATION_PARAMETERS, is_admissible +from rubem.calibration.runner import EVALUATION_COLUMNS, CalibrationSettings, calibrate +from rubem.configuration.model_configuration import ModelConfiguration + +# SciPy is the optional ``rubem[calibration]`` extra; an environment without it +# collects this file and skips it instead of failing to import. +pytest.importorskip("scipy.optimize") + +pytestmark = [pytest.mark.dataset, pytest.mark.slow] + +BASIN = "ipojuca_basin" +"""The basin this test calibrates: its manifest is ``.json``.""" + +PATH_SEGMENT = f"{BASIN}/" +"""The segment every input path of the manifest carries. + +The manifests are written by the machine that extracted the datasets and hold +absolute paths of that machine, so every path is rebased onto the local copy by +replacing whatever precedes this segment. +""" + +STEPS = 6 +"""Monthly steps of the window the test simulates, from 01/01/2000 to 01/06/2000.""" + +# Five admissible members around the parameter set of the manifest, the second +# one far from it. SciPy replaces the first member with x0, so the configuration +# under calibration is always part of the initial population, and the search has +# one generation to find it again. +INIT = np.array( + [ + [4.415, 0.078, 0.51, 0.12, 5.375, 0.581, 0.922, 0.307], + [9.500, 0.900, 0.05, 0.85, 9.000, 0.950, 0.150, 0.900], + [3.500, 0.150, 0.45, 0.20, 4.500, 0.500, 0.800, 0.250], + [5.500, 0.050, 0.55, 0.10, 6.500, 0.650, 0.950, 0.400], + [4.000, 0.200, 0.40, 0.30, 5.000, 0.450, 0.700, 0.350], + ] +) + + +def _rebase(value, root: Path): + """Return the manifest with every input path anchored on the local copy. + + :param value: A manifest document, one of its sections or one of its values. + :param root: The directory the extracted datasets live in. + :type root: pathlib.Path + """ + if isinstance(value, dict): + return {key: _rebase(item, root) for key, item in value.items()} + if isinstance(value, str) and PATH_SEGMENT in value: + return str(root / value[value.index(PATH_SEGMENT) :]) + return value + + +def _input_paths(value, root: Path): + """Yield every path of a rebased document that lies under ``root``.""" + if isinstance(value, dict): + for item in value.values(): + yield from _input_paths(item, root) + elif isinstance(value, str) and value.startswith(str(root)): + yield value + + +def _read_evaluations(path: Path) -> list[dict[str, str]]: + """Return the rows of an ``evaluations.csv`` table, header first. + + A local copy of the reader of ``tests/unit/calibration/test_runner.py``: + importing one test module from another would run the module-level + ``importorskip`` of the other one during the collection of this one, for + three lines of CSV reading. + """ + with path.open(encoding="utf-8", newline="") as table: + return list(csv.DictReader(table)) + + +def _fixed_ldd(document: dict, destination: Path) -> str: + """Generate the drainage network of the basin once and return its path. + + The network is derived from the digital elevation model of the basin with + the parameters the model itself uses when no network is configured, and is + written so that every run of the calibration routes the runoff along the + very same network. + """ + import pcraster as pcr + + pcr.setclone(document["RASTERS"]["clone"]) + dem = pcr.readmap(document["RASTERS"]["dem"]) + pcr.report(pcr.lddcreate(dem, 1e31, 1e31, 1e31, 1e31), str(destination)) + return str(destination) + + +def test_the_calibration_recovers_the_parameters_of_the_ipojuca_basin( + dataset_dir, tmp_path, monkeypatch +): + manifest = dataset_dir / f"{BASIN}.json" + if not manifest.is_file(): + pytest.skip(f"{manifest} is not there") + + # The inputs are read where they were extracted, so the persistent auxiliary + # metadata of GDAL would annotate every raster of the dataset with a sidecar + # file of its statistics, next to the raster itself. The test writes nothing + # outside its own temporary directory; the workers inherit the variable. + monkeypatch.setenv("GDAL_PAM_ENABLED", "NO") + + assert [tuple(row) for row in INIT if not is_admissible(row)] == [], ( + "every member of the initial population must be admissible, or the search " + "would spend its generation on candidates it never runs" + ) + + document = _rebase(json.loads(manifest.read_text(encoding="utf-8")), dataset_dir) + missing = [path for path in _input_paths(document, dataset_dir) if not Path(path).exists()] + assert missing == [], "the manifest was rebased onto the local copy of the datasets" + + output_dir = tmp_path / "output" + output_dir.mkdir() + document["DIRECTORIES"]["output"] = str(output_dir) + document["SIM_TIME"] = {"start": "01/01/2000", "end": "01/06/2000"} + document["RASTERS"]["ldd"] = _fixed_ldd(document, tmp_path / "ldd.map") + + config_file = tmp_path / f"{BASIN}.json" + config_file.write_text(json.dumps(document, indent=2), encoding="utf-8") + parameters = ModelConfiguration(config_file, validate_input=False).calibration_parameters + expected = parameters.model_dump() + + written = Model.from_file(config_file).run().time_series["arn"][0] + observed = tmp_path / "observed.csv" + shutil.copyfile(written, observed) + assert len(observed.read_text(encoding="utf-8").splitlines()) == STEPS + 1, ( + "the header and one row per monthly step of the window" + ) + + temp_dir = tmp_path / "tmp" + result = calibrate( + config_file, + observed, + tmp_path / "run", + CalibrationSettings( + init=INIT, + maxiter=1, + popsize=5, + workers=2, + seed=1, + temp_dir=str(temp_dir), + ), + ) + + # The observed series is what the configuration itself wrote, so the global + # optimum of the search is that configuration, with an efficiency of 1. + assert result.best_objective == 0.0 + assert result.best_nse == pytest.approx(1.0, abs=1e-12) + assert result.best_parameters == pytest.approx(expected, abs=1e-12) + assert set(result.best_parameters) == set(CALIBRATION_PARAMETERS) + + rows = _read_evaluations(result.evaluations_csv) + assert list(rows[0]) == list(EVALUATION_COLUMNS) + assert len(rows) == result.evaluations + assert [row["error"] for row in rows if row["error"]] == [] + + summary = json.loads(result.result_json.read_text(encoding="utf-8")) + assert summary["nfev"] == result.evaluations + assert summary["population_size"] == len(INIT) + + assert result.calibrated_config.name == f"{BASIN}-calibrated.json" + calibrated = ModelConfiguration(result.calibrated_config, validate_input=False) + assert calibrated.calibration_parameters.model_dump() == pytest.approx( + result.best_parameters, abs=1e-12 + ) + + assert multiprocessing.active_children() == [] + assert list(temp_dir.iterdir()) == [] From 9914a7472e122725ff72811871dca87a8528ced6 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Wed, 16 Sep 2026 18:30:13 -0300 Subject: [PATCH 10/12] Document the calibration method Add the calibration page: what the command fits, the objective and its masks and guards, the free parameters with their bounds and the derived weight, the differential evolution settings and the evaluation budget, the observed series formats and the fixed drainage network a calibration on the accumulated runoff needs, the artifacts, the process model and hints for large machines; link it from the user guide and the API page. --- doc/source/api.rst | 31 +++ doc/source/calibration.rst | 547 +++++++++++++++++++++++++++++++++++++ doc/source/index.rst | 1 + doc/source/userguide.rst | 5 +- 4 files changed, 582 insertions(+), 2 deletions(-) create mode 100644 doc/source/calibration.rst diff --git a/doc/source/api.rst b/doc/source/api.rst index 96344b8a..d432552f 100644 --- a/doc/source/api.rst +++ b/doc/source/api.rst @@ -122,6 +122,37 @@ Importing ``rubem.api`` does not require PCRaster or GDAL; running the model does. Without them, loading a configuration or running raises ``ImportError`` with the same installation guidance the command line prints. +Calibration +----------- + +``rubem calibrate`` is built on this API: each of its workers calls +``Model.run()`` in its own process, on a configuration derived from the one +being calibrated. The same search is available from Python as +``rubem.calibration.runner.calibrate``, which takes the configuration file, the +observed series, the run directory and a ``CalibrationSettings``: + +.. code-block:: python + + from rubem.calibration.runner import CalibrationSettings, calibrate + + if __name__ == "__main__": + result = calibrate( + "config.json", + "observed.csv", + "calibration", + CalibrationSettings(seed=42, maxiter=20, popsize=5), + ) + print(result.best_parameters, result.best_nse) + +The entry point guard is required: the calibration starts its workers with the +``spawn`` method, which imports the main module of the caller in every one of +them. SciPy is needed as well, from the ``rubem[calibration]`` extra. + +``rubem.calibration`` is **not** part of the stable surface yet: only +``rubem.api`` is, and the calibration package may change without a minor +version bump until it is promoted. The method behind the command is documented +in :doc:`Calibration `. + Logging ------- diff --git a/doc/source/calibration.rst b/doc/source/calibration.rst new file mode 100644 index 00000000..ea0643d9 --- /dev/null +++ b/doc/source/calibration.rst @@ -0,0 +1,547 @@ +Calibration +=========== + +.. role:: raw-html(raw) + :format: html + +The ``rubem calibrate`` command fits the calibration parameters of the model to +an observed streamflow series with a differential evolution search +[STORN1997]_. Every candidate the search proposes is a complete simulation: +the command builds the configuration of the candidate, runs it, samples the +resulting time series at the stations and scores it against the observations. + +The options of the command are listed in the user guide, under +:ref:`userguide:Running RUBEM`. This page documents what happens between them: +the objective, the parameters and their bounds, the settings of the search and +how to budget a run, the inputs it needs, what it writes, and how it uses the +machine it runs on. + +What the calibration fits +------------------------- + +The model has nine calibration parameters, the ``CALIBRATION`` section of the +configuration file: :math:`\alpha`, :math:`b`, :math:`w_1`, :math:`w_2`, +:math:`w_3`, :math:`RCD`, :math:`f`, :math:`\alpha_{GW}` and :math:`x`. They +are described in :ref:`overview:Calibration and Validation`. The search covers +eight of them; the slope factor weight :math:`w_3` is not searched but derived +from the other two weights, because the three weights of the potential runoff +coefficient must add up to 1. + +The observations are a streamflow series at the sample stations of the +configuration, one column per station. By default the simulated counterpart is +``arn``, the accumulated total runoff, in +:raw-html:`m3s-1`, which is the variable an observed +streamflow record is directly comparable with. ``--variable`` selects another +output variable, one of ``itp``, ``bfw``, ``srn``, ``eta``, ``lfw``, ``rec``, +``smc``, ``rnf`` and ``arn``; note that ``rnf`` is the total runoff of a single +cell, in millimeters, not a discharge, so an observed series compared with it +must be in the same quantity. + +Each evaluation runs the whole simulation period of the configuration. The +calibration never changes the period, the inputs or the aggregation: of what +the results depend on, it changes the nine parameters and nothing else. + +The objective function +---------------------- + +The efficiency +`````````````` + +The agreement between a simulated and an observed series at one station is the +Nash-Sutcliffe efficiency [MCCUEN2006]_: + +.. math:: + :label: calibration-nse + :nowrap: + + \[NSE = 1 - \frac{\sum_{t}{\left(Q_{sim,t} - Q_{obs,t}\right)^2}}{\sum_{t}{\left(Q_{obs,t} - \overline{Q_{obs}}\right)^2}}\] + +where: + +- :math:`Q_{sim,t}` – simulated value of the station at time step :math:`t`; +- :math:`Q_{obs,t}` – observed value of the station at time step :math:`t`; +- :math:`\overline{Q_{obs}}` – mean of the observed values used in the sum. + +Both sums run over the same pairs, and the mean is the mean of the observed +values of those pairs: a step dropped because the simulation has no value for +it is dropped from the mean as well, so numerator and denominator always +describe the same sample. + +Which values count +`````````````````` + +A pair counts only when both of its members count. A value does not count when +it is not finite (``NaN`` and the infinities), when it is ``-9999``, the gap +marker of the model's own tables and of the usual observation tables, or when +it is at or above ``1e30``. + +That last threshold is how the PCRaster missing value ``1e31`` is masked. The +value does not round-trip exactly through ``Float32`` — written and read back +it becomes ``9.999999...e30`` — so the mask rejects everything at or above +``1e30`` instead of comparing with the constant. No streamflow or runoff of a +real basin comes anywhere near that magnitude. An empty cell, or a cell that is +not a number at all, is read as ``NaN`` and masked the same way. + +Stations and the average +```````````````````````` + +The two series are aligned on the station ids they share and on the time steps +they share; a station or a step present on one side only is ignored. + +A station has no efficiency when fewer than two of its pairs count, or when all +of its valid observations are equal: a constant record has no variability +for the efficiency to explain and would divide by zero. Such a station is +reported as having no efficiency and is left out of the average, instead of +dragging it. The efficiency of a candidate is the mean over the stations that +do have one. + +A mismatch between the observations and the configuration — no station id in +common, no time step in common after the spin-up, or no shared station with an +efficiency — is not a poor candidate. It is not caught before the search +either: the comparison happens inside the evaluation, after that evaluation has +run the whole simulation. The evaluation is then recorded with the message of +the error in its ``error`` column and with the objective ``1e30``, like any +other failure, and the search carries on. The calibration ends with an error +only once no evaluation of the whole run has succeeded. + +.. warning:: + + A wrong station id is therefore paid for in model runs before it is + reported. The per-generation log line is the early sign: a best objective of + ``1e+30`` from the first generation on means that every candidate failed. + Check the observed file against a run of the configuration first, over a few + time steps, as `Running on a cluster or a large machine`_ describes. + +The spin-up window +`````````````````` + +``--spinup-steps`` excludes the leading time steps of the simulation from the +comparison: every step whose number is at most ``--spinup-steps`` is dropped +from both series before they are aligned. The simulation still runs from its +first step — the point is to let the storages of the model fill before their +effect is scored, not to shorten the run. + +From the efficiency to the objective +```````````````````````````````````` + +The differential evolution minimizes, so the efficiency is turned into a cost: + +.. math:: + :label: calibration-objective + :nowrap: + + \[FO = 1000 \cdot \left(100 \cdot \left(1 - NSE\right)\right)^2\] + +A perfect simulation, :math:`NSE = 1`, gives :math:`FO = 0`; a simulation no +better than the mean of the observations, :math:`NSE = 0`, gives +:math:`FO = 10^7`; and the cost grows quadratically as the efficiency falls +further. + +A candidate that is never run — one the admissibility check rejects — and a +candidate whose run or evaluation fails are both recorded with the objective +``1e30``. The worst value a real run can produce is far below it (an efficiency +of :math:`-100`, already an absurd simulation, gives about +:math:`1.0 \times 10^{11}`), so a rejected or failed candidate always ranks +behind every candidate that was actually evaluated. The value is finite on +purpose: it sorts, it is written to JSON, and it appears in the table of +evaluations like any other objective. + +Free parameters and their bounds +-------------------------------- + +The eight searched parameters, in the order of the decision vector, with the +bounds the application settings declare (:file:`rubem/appsettings.json`, +section ``value_ranges.variables``). These are the same ranges the +configuration loader validates against, so any candidate inside the bounds is a +configuration the model accepts. + +.. list-table:: Searched parameters + :header-rows: 1 + :widths: 20 20 40 20 + + * - Configuration key + - Symbol + - Description + - Bounds + * - ``alpha`` + - :math:`\alpha` + - Interception parameter + - :math:`0.01 \leq \alpha \leq 10` + * - ``b`` + - :math:`b` + - Rainfall intensity coefficient + - :math:`0.01 \leq b \leq 1` + * - ``w_1`` + - :math:`w_1` + - Land use factor weight + - :math:`0 \leq w_1 \leq 1` + * - ``w_2`` + - :math:`w_2` + - Soil factor weight + - :math:`0 \leq w_2 \leq 1` + * - ``rcd`` + - :math:`RCD` + - Regional consecutive dryness level + - :math:`1 \leq RCD \leq 10` + * - ``f`` + - :math:`f` + - Flow direction factor + - :math:`0.01 \leq f \leq 1` + * - ``alpha_gw`` + - :math:`\alpha_{GW}` + - Baseflow recession coefficient + - :math:`0.01 \leq \alpha_{GW} \leq 1` + * - ``x`` + - :math:`x` + - Flow recession coefficient + - :math:`0 \leq x \leq 1` + +.. note:: + + The configuration file writes the rainfall intensity coefficient as ``b``, + and so does the calibrated configuration the run writes; + :file:`evaluations.csv` and :file:`result.json` spell it ``beta``, which is + the name the parameter carries inside the package. They are the same + parameter. + +The ninth parameter is derived rather than searched: + +.. math:: + :label: calibration-w3 + :nowrap: + + \[w_3 = 1 - \left(w_1 + w_2\right)\] + +so the search is eight-dimensional under the linear constraint +:math:`w_1 + w_2 \leq 1`. The constraint is handed to the optimizer, which +checks it before the objective and never spends a model run on a candidate +whose slope factor weight would come out negative. The derived weight must also +lie inside its own range, :math:`0 \leq w_3 \leq 1`; a candidate for which it +does not is rejected without a run and recorded with the error +``inadmissible``. Both checks are guards rather than a part of the search: the +optimizer keeps its trial vectors inside the bounds and drops the ones that +violate the constraint, so it does not propose an inadmissible candidate of its +own. + +The starting point of the search, ``x0``, is the parameter set of the +configuration being calibrated. The optimizer places it in the initial +population, so the configuration's own parameters are always among the +candidates of the first generation and a calibration can never return something +worse than what it started from without the table showing why. + +The search +---------- + +Settings +```````` + +The search is SciPy's ``differential_evolution`` with the following settings, +which the command does not expose because they describe the method rather than +the run: + +``strategy="best1exp"`` + The trial vector is built from the best member of the population plus a + scaled difference of two random members, with exponential crossover. A + best-based strategy exploits the best candidate found so far and converges + in fewer generations than the random-base ones [STORN1997]_, at some cost in + exploration, which the dithered mutation below offsets. Fewer generations is + what matters when one evaluation is a full simulation. + +``init="sobol"`` + The initial population is a Sobol' sequence over the bounds instead of a + random or Latin hypercube sample: it covers the eight-dimensional box far + more evenly, so the initial sample is not wasted on clusters and gaps. + +``mutation=(0.5, 1.0)`` + A pair, not a number: the mutation constant is drawn anew for every + generation between the two values (*dithering*). It widens the range of step + sizes the search takes without any tuning and reduces the chance of + stagnating in a local optimum. + +``recombination=0.7`` + The probability with which a trial vector keeps a component of the mutant + rather than of its parent. A high value suits parameters that interact, as + the three weights and the runoff coefficients do. + +``polish=False`` + No local refinement of the best candidate at the end. The local search would + spend model runs estimating gradients of an objective that is read back from + a simulated table, which does not provide them reliably. + +``updating="deferred"`` + The whole generation is evaluated before the population is updated. This is + what makes the generation a batch that can be spread over the worker + processes; the immediate alternative is inherently serial. + +``--seed`` seeds the random generator of the search: the same seed gives the +same initial population and the same random draws of the mutation and of the +crossover. That alone does not make two calibrations comparable candidate by +candidate, because from the first generation on the selection reads the +objectives, so the candidates a run proposes depend on the values the +simulations returned. Two runs follow the same path only when the simulation +itself is reproducible — see `A fixed drainage network`_. + +The evaluation budget +````````````````````` + +The number of population members is not ``--popsize``: it is + +.. math:: + :label: calibration-members + :nowrap: + + \[members = \max\left(5,\ popsize \times 8\right)\] + +with 8 the number of free parameters, rounded up to the next power of two, +because a Sobol' sequence is balanced only over a power-of-two sample. Every +generation evaluates one candidate per member, and the initial population is +one more round of evaluations on top of ``--maxiter`` generations, so + +.. math:: + :label: calibration-runs + :nowrap: + + \[runs \leq members \times \left(maxiter + 1\right)\] + +With the defaults, ``--popsize 15`` gives :math:`15 \times 8 = 120`, rounded up +to **128 members**, and ``--maxiter 100`` gives at most +:math:`128 \times 101 = 12\,928` model runs. A smaller search, ``--popsize 5``, +gives :math:`5 \times 8 = 40`, rounded up to **64 members**, so +``--popsize 5 --maxiter 20`` costs at most :math:`64 \times 21 = 1\,344` runs. + +Two things make the real count lower. Members whose weights violate +:math:`w_1 + w_2 \leq 1` are skipped by the optimizer without a model run and +leave no record. And the search stops as soon as the objective values of the +population agree: SciPy compares their standard deviation with its default +tolerance of 1% of their mean, and skips the comparison while any member of +the population is still infeasible. A search that converges this way stops +before ``--maxiter`` is reached. + +The number of members and the resulting upper bound are logged when the +calibration starts, so the budget can be checked before the machine is +committed to it: + +.. code-block:: console + + Calibrating 8 free parameter(s) with 128 population member(s) per generation and + at most 100 generation(s): up to 12928 model run(s), on 15 worker process(es). + +The size of the population the search actually built is recorded as +``population_size`` in :file:`result.json`, together with the number of +evaluations and of generations it really spent. + +Inputs +------ + +The observed series +``````````````````` + +``--observed`` takes one file, in either of two layouts, recognized from its +content. + +The **CSV layout** is the one the model writes for its own time series tables: +``;`` as the separator, a header whose first cell is ``0`` and whose remaining +cells are the station ids, and one row per time step carrying the step number +in the first column and one value per station: + +.. code-block:: console + + 0;1;2;3 + 1;12.4;3.1;0.8 + 2;15.9;4.0;1.2 + 3;-9999;4.4;1.0 + +The **PCRaster TSS layout** is accepted as well, both in its documented form (a +title line, the number of columns, one line per column name, then the rows) and +in the header-less form the model itself writes, whose stations are then +numbered ``1``, ``2``, ... in column order. + +.. note:: + + The station ids must be the ids of the sample locations of the + configuration. The easiest way to get the layout and the ids right is to run + the configuration once and use the table the run wrote as the template of + the observed file: :file:`tss_arn.csv` for the default ``point`` + aggregation, :file:`tss_arn_subcatchment.csv` or :file:`tss_arn_zones.csv` + for the other two. With ``zones`` the ids are the renumbered columns + ``1..N``, and the correspondence with the original zone ids is written to + :file:`zones_mapping.csv` in the output directory. + +The configuration must name the raster the aggregation samples from — the +sample locations raster for ``point`` and ``subcatchment``, the zones raster +for ``zones``. A configuration that does not name the raster its aggregation +needs is refused before any run: there would be no station to compare. + +A fixed drainage network +```````````````````````` + +When the configuration does not provide a Local Drain Direction raster, the +model derives one from the digital elevation model with PCRaster's +``lddcreate``. On the real basin DEMs that operation is not deterministic: two +runs of the same configuration produce drainage networks that differ in a few +cells, and the accumulated runoff ``arn``, which is routed over the network, +therefore differs between them. The other eight output variables are computed +cell by cell and are unaffected. + +A calibration on ``arn`` without a fixed network would compare series routed +over different networks from one evaluation to the next, and the differences +between candidates would be partly noise. **Calibrating on** ``arn`` +**requires** ``RASTERS.ldd`` **to name a fixed LDD raster.** Generate it once, +with the DEM of the basin, and keep it with the other inputs: + +.. code-block:: python + + import pcraster as pcr + + pcr.setclone("/basins/ipojuca/input/maps/clone.map") + dem = pcr.readmap("/basins/ipojuca/input/maps/dem.map") + ldd = pcr.lddcreate(dem, 1e31, 1e31, 1e31, 1e31) + pcr.report(ldd, "/basins/ipojuca/input/maps/ldd.map") + +.. code-block:: json + + { + "RASTERS": { + "ldd": "/basins/ipojuca/input/maps/ldd.map", + }, + } + +The same applies beyond ``arn`` when the time series are aggregated over +subcatchments: the subcatchment of each sample location is delineated over the +drainage network, so the areas the series are averaged over move with it, for +every variable. + +Artifacts of a run +------------------ + +The three artifacts of the run directory, ``-o``/``--run-dir``: + +:file:`evaluations.csv` + One row per evaluation, with the columns ``id``, ``pid``, ``alpha``, + ``beta``, ``w_1``, ``w_2``, ``w_3``, ``rcd``, ``f``, ``alpha_gw``, ``x``, + ``nse``, ``objective``, ``elapsed_seconds`` and ``error``. ``id`` identifies + the evaluation and ``pid`` the process that made it; ``nse`` is empty for a + candidate that produced no efficiency; ``error`` is empty for a successful + evaluation, ``inadmissible`` for a candidate rejected without a run, and the + type and message of the exception for a run that failed. The rows are + ordered by the candidate they evaluated, not by the moment they were + written, so that the table does not depend on how the parallel workers + happened to finish. + +:file:`result.json` + The summary: ``best_parameters`` (the nine parameters, the derived ``w_3`` + included), ``best_nse``, ``best_objective``, ``nfev`` and ``nit`` (the + evaluations and generations the search spent), ``success`` and ``message`` + from the optimizer, ``population_size``, and a ``settings`` object with + ``variable``, ``spinup_steps``, ``maxiter``, ``popsize``, ``seed``, + ``workers``, ``temp_dir``, ``init``, ``strategy``, ``mutation``, + ``recombination`` and ``polish``. The number of workers and the temporary + directory are the values the run resolved and used, not the defaults left + unset on the command line, so the summary describes a calibration that can + be repeated. + +:file:`-calibrated.json` + The calibrated configuration: the input document with the best parameters + in place, written in the format the input was written in, legacy or + format 1.0. Its paths are the ones the loader resolved, which are absolute + even when the input file wrote them relative to its own directory; moving + the file to another machine therefore means fixing the paths. + +Next to them, the :file:`evaluations/` directory holds one JSON record per +evaluation, written by the worker that made it. A record carries what the CSV +row carries and, in addition, ``station_nse``, the efficiency of each +individual station — which is where to look when the mean hides one station +behaving differently from the others. + +.. warning:: + + A run directory that already holds evaluation records of an earlier + calibration is refused before the search starts: consolidating two runs into + one table would mix them. Use an empty run directory for every calibration, + and keep the finished ones. + +Process model of the calibration +-------------------------------- + +The parent process loads the configuration once and validates its input files +once. The workers do not validate them again: they do not change during the +calibration, and re-reading every raster of the basin per evaluation would +dominate the cost. + +Every evaluation then runs in a **fresh interpreter**. PCRaster keeps its state +process-wide — ``setclone`` is global to the process and the raster memory of a +run is not reclaimed — so a worker reused for a second evaluation would carry +the state of the first. The pool is therefore started with the ``spawn`` method +and configured for one task per worker: a worker starts, evaluates one +candidate and exits. The generation is spread over as many such workers as +``--workers`` allows, which is also what makes the calibration parallel at all. + +Each worker builds the configuration of its candidate from the document of the +calibrated configuration: the nine parameters become the candidate's, the +output directory becomes a temporary directory of its own under ``--temp-dir``, +every raster series is disabled (an evaluation reads no raster of a previous +run, and writing them would dominate the cost of the run) and the time series +are reduced to the calibrated variable, as CSV. The aggregation the user +configured is kept, since it decides which areas the observed stations +correspond to. The temporary directory is removed when the evaluation ends, +whether it succeeded or not. + +The workers are silenced: a null handler on the root logger keeps the warnings +of a configuration that is not validated again out of the terminal, and +standard output is pointed at the null device so that the newline the PCRaster +framework writes when an interpreter exits is not printed once per evaluation. +Nothing is lost: a failure of an evaluation is recorded in that evaluation's +JSON record and in its row of the table, and the parent logs one summary at the +end saying how many evaluations failed and what the first failure was. One +failing candidate never ends a calibration; it is ranked behind every candidate +that was evaluated. + +Running on a cluster or a large machine +--------------------------------------- + +The cost of a calibration is the number of model runs times the cost of one +run, divided by the number of workers. Measure the second factor before +committing a machine: one ``rubem run`` of the configuration is an upper bound +on one evaluation, which writes no raster series and does not validate the +inputs again, and so costs a little less. On the Ipojuca basin, 228 monthly +time steps, an evaluation is of the order of half a minute. At 25 seconds +each, a default search on that basin, up to 12,928 runs on 15 workers, is on +the order of six hours; ``--popsize 5 --maxiter 20``, up +to 1,344 runs, is well under an hour. + +``--workers`` + Defaults to one less than the number of CPUs, which leaves a core for the + parent. On a shared machine, or under a batch scheduler, set it to the + number of cores the job was actually given rather than letting it read the + whole machine. Each worker holds one model run in memory, plus the parent, + so the memory of the job is roughly ``--workers`` times the memory of a + single ``rubem run`` of the same configuration; that, and not the core + count, is usually what limits the number of workers on a large basin. + +``--temp-dir`` + The per-evaluation output directories are created here. Point it at fast + local storage of the node — the scratch directory of the job, not a shared + network filesystem: every evaluation creates a directory, writes a small + table into it and removes it again, and a network filesystem turns that into + the slowest part of the evaluation. The run directory itself may stay on + shared storage; it receives one small JSON record per evaluation. + +Shorten the simulation period first + An exploratory calibration over a few months of ``SIM_TIME``, with a small + ``--popsize`` and ``--maxiter``, costs a fraction of the full search and + still shows whether the observed series, the station ids and the spin-up are + set up as intended. Move to the full period only once a short run finishes + and produces a sensible efficiency. + +Continue from a calibrated configuration + :file:`-calibrated.json` is an ordinary configuration file. Passing + it back to ``rubem calibrate`` as ``-c`` starts a second search from the + parameters the first one found, since the starting point of a search is the + parameter set of the configuration it is given — useful to refine a result + with a different seed, a longer period or a larger population. + +Keep the run directory + :file:`evaluations.csv` and :file:`result.json` are the record of how a + parameter set was obtained: the settings, the seed, the number of + evaluations and the whole response surface the search sampled. They are + small, and they are what makes a published calibration reproducible. diff --git a/doc/source/index.rst b/doc/source/index.rst index d2eff553..36615ffc 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -32,6 +32,7 @@ All development activity is coordinated via the `RUBEM Github page -calibrated.json`` into the ``--run-dir``. The method itself (the -objective, the parameter bounds, the derivation of ``w3`` and how to budget the -number of evaluations) is documented separately. +objective, the parameter bounds, the derivation of ``w3``, how to budget the +number of evaluations and what a run needs from the configuration) is +documented in :doc:`Calibration `. .. code-block:: console From acea648586fdb1b1c6b3c9749e86051ac017ad34 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Wed, 16 Sep 2026 18:30:14 -0300 Subject: [PATCH 11/12] Record the calibration documentation and the dataset marker in the changelog --- doc/source/changelog.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 385dccf9..da34c43b 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -13,6 +13,14 @@ The format follows `Keep a Changelog `__. Added ````` +- Added the calibration page of the documentation (the objective, the free + parameters and their bounds, the derivation of ``w3``, the differential + evolution settings and the evaluation budget, the observed series formats, + the fixed drainage network a calibration on ``arn`` needs, the artifacts, the + process model and hints for large machines), the ``dataset`` pytest marker + with the ``RUBEM_DATASET_DIR`` convention for tests that need the published + datasets, and a calibration smoke test on the Ipojuca basin + (`#345 `__). - Added ``rubem calibrate``: a differential evolution (SciPy) over the eight free calibration parameters, ``w3`` derived from the other two weights, minimizing ``1000 (100 (1 - NSE))^2`` on the station series of one output variable From 5d0d77f251938f704b9b7598d4ae7559acb74afa Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Wed, 16 Sep 2026 18:51:20 -0300 Subject: [PATCH 12/12] Describe the refusal of an incomparable observed series before the search The calibration page said every mismatch between the observations and the configuration surfaced inside the evaluations; the runner now refuses a series with no shared time step or no shared station before the search starts, and only a series whose shared stations all lack an efficiency still fails each evaluation after its simulation. --- doc/source/calibration.rst | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/doc/source/calibration.rst b/doc/source/calibration.rst index ea0643d9..a8bb24cd 100644 --- a/doc/source/calibration.rst +++ b/doc/source/calibration.rst @@ -95,22 +95,30 @@ reported as having no efficiency and is left out of the average, instead of dragging it. The efficiency of a candidate is the mean over the stations that do have one. -A mismatch between the observations and the configuration — no station id in -common, no time step in common after the spin-up, or no shared station with an -efficiency — is not a poor candidate. It is not caught before the search -either: the comparison happens inside the evaluation, after that evaluation has -run the whole simulation. The evaluation is then recorded with the message of -the error in its ``error`` column and with the objective ``1e30``, like any -other failure, and the search carries on. The calibration ends with an error -only once no evaluation of the whole run has succeeded. +A mismatch between the observations and the configuration is not a poor +candidate, and the two kinds of it that can be told before any model run are +refused before the search starts: an observed series with no time step in +common with the simulated steps that survive the spin-up, and one whose station +ids match none of the ids of the sample locations raster (an observed station +the configuration does not sample is reported once and ignored). With the +``zones`` aggregation the station ids are only known once a run has written +``zones_mapping.csv``, so that check is skipped. + +What cannot be told in advance surfaces inside the evaluations: a series whose +shared stations all lack an efficiency, because every observation is missing or +constant, fails each evaluation after its simulation has run. The evaluation is +then recorded with the message of the error in its ``error`` column and with +the objective ``1e30``, like any other failure, and the search carries on; the +calibration ends with an error only once no evaluation of the whole run has +succeeded. .. warning:: - A wrong station id is therefore paid for in model runs before it is - reported. The per-generation log line is the early sign: a best objective of - ``1e+30`` from the first generation on means that every candidate failed. - Check the observed file against a run of the configuration first, over a few - time steps, as `Running on a cluster or a large machine`_ describes. + Such a failure is paid for in model runs before it is reported. The + per-generation log line is the early sign: a best objective of ``1e+30`` from + the first generation on means that every candidate failed. Check the observed + file against a run of the configuration first, over a few time steps, as + `Running on a cluster or a large machine`_ describes. The spin-up window ``````````````````