From 0a3e08ff8f52be69e3d9797495a0354016ab9b8a Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Fri, 18 Sep 2026 16:12:41 -0300 Subject: [PATCH 1/9] 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, and give rubem._deps the guard of the optional dependency, whose guidance names the extra. 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 +++++- rubem/_deps.py | 50 +++++++++++++++++++++++++++++++++++++++ tests/unit/test_deps.py | 51 ++++++++++++++++++++++++++++++++++++++++ 6 files changed, 144 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", 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/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 5d72e938ca21d6949a6d0ab9322fbb4425ba8853 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Fri, 18 Sep 2026 16:12:41 -0300 Subject: [PATCH 2/9] Add the calibration objective, the station metrics and the series readers rubem.calibration.objective computes the Nash-Sutcliffe efficiency over the pairs valid in both series (NaN, negative values, -9999 and the PCRaster missing value are gaps; fewer than two pairs or a constant observed record give no value), the objective 1000 (100 (1 - NSE))^2 and, per shared station, the pairs, the means and standard deviations of both series, the Pearson correlation, the RMSE and the NSE. evaluate_series aligns a simulated and an observed series on their shared time steps and station ids after a spin-up window and averages the NSE over a chosen selection of stations while measuring the others. read_series reads the table layout of the model's time series ("0;;...") and PCRaster time series with their header; a file without a header is refused naming the two accepted layouts, and a byte order mark is ignored. --- rubem/calibration/__init__.py | 12 + rubem/calibration/objective.py | 627 ++++++++++++++++++++ tests/unit/calibration/__init__.py | 0 tests/unit/calibration/test_objective.py | 695 +++++++++++++++++++++++ 4 files changed, 1334 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..2272fb49 --- /dev/null +++ b/rubem/calibration/objective.py @@ -0,0 +1,627 @@ +"""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 goodness of fit of each +station and turns the mean Nash-Sutcliffe efficiency (NSE) into the objective +the optimizer minimizes. + +Nothing here imports SciPy: the efficiency, the remaining statistics and the +objective are evaluated with NumPy alone. +""" + +import csv +import logging +from collections.abc import Sequence +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. + +``-9999`` is the one the model writes, and it is negative, so the sign rule of +:func:`valid_mask` would drop it on its own; it is named here because a table +that writes it is saying "gap" and not "an implausible discharge". +""" + +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 negative, is not one of + :data:`MISSING_VALUES` and is below ``1e30``, which covers the PCRaster + missing value in any of its round-trips. + + Every quantity the two series carry is a flux or a storage and cannot be + negative, so a negative entry is not a measurement: the observation tables + of the sample stations write the gaps of a gauge as ``-9999``, as ``-1`` or + as another negative marker of their own, and the mask drops all of them + alike. Zero is a value, not a gap: a station may well record no flow. + + :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) & (array >= 0.0) + 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 + + +@dataclass(frozen=True) +class StationMetrics: + """The goodness of fit of one station over the pairs both series have. + + Every field but :attr:`pairs` is ``None`` where the statistic is not + defined: a mean and a root mean squared error need one pair, a standard + deviation with ``ddof=1`` and an efficiency need two, and a correlation and + an efficiency need the series they divide by to vary. + + :param pairs: Number of time steps at which both series carry a value the + mask of :func:`valid_mask` accepts. + :type pairs: int + + :param mean_observed: Mean of the observed values of those pairs. + :type mean_observed: float | None + + :param std_observed: Standard deviation of the observed values of those + pairs, with one degree of freedom (``ddof=1``). + :type std_observed: float | None + + :param mean_simulated: Mean of the simulated values of those pairs. + :type mean_simulated: float | None + + :param std_simulated: Standard deviation of the simulated values of those + pairs, with one degree of freedom (``ddof=1``). + :type std_simulated: float | None + + :param r: Pearson correlation coefficient of the two series over those + pairs, ``None`` when either of them is constant. + :type r: float | None + + :param rmse: Root mean squared error over those pairs. + :type rmse: float | None + + :param nse: Nash-Sutcliffe efficiency, as :func:`nash_sutcliffe` defines it. + :type nse: float | None + """ + + pairs: int + mean_observed: float | None + std_observed: float | None + mean_simulated: float | None + std_simulated: float | None + r: float | None + rmse: float | None + nse: float | None + + +def station_metrics(simulated: np.ndarray, observed: np.ndarray) -> StationMetrics: + """Return the goodness of fit of one station. + + Every statistic is computed over the pairs whose two members pass + :func:`valid_mask`, the same sample :func:`nash_sutcliffe` uses, so the + efficiency and the statistics that explain it always describe one 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 statistics of the station. + :rtype: StationMetrics + + :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) + pairs = int(np.count_nonzero(mask)) + if pairs == 0: + return StationMetrics(0, None, None, None, None, None, None, None) + + simulated_valid = simulated[mask] + observed_valid = observed[mask] + mean_simulated = float(np.mean(simulated_valid)) + mean_observed = float(np.mean(observed_valid)) + rmse = float(np.sqrt(np.mean((simulated_valid - observed_valid) ** 2))) + if pairs < 2: + return StationMetrics(pairs, mean_observed, None, mean_simulated, None, None, rmse, None) + + return StationMetrics( + pairs=pairs, + mean_observed=mean_observed, + std_observed=float(np.std(observed_valid, ddof=1)), + mean_simulated=mean_simulated, + std_simulated=float(np.std(simulated_valid, ddof=1)), + r=_correlation(simulated_valid, observed_valid), + rmse=rmse, + nse=nash_sutcliffe(simulated_valid, observed_valid), + ) + + +def _correlation(simulated: np.ndarray, observed: np.ndarray) -> float | None: + """Return the Pearson correlation of two samples, or ``None``. + + The coefficient is undefined when either sample is constant, which leaves + its own spread at zero and the quotient without a denominator. Constancy is + tested on the values themselves, since the sum of the squared deviations of + a constant sample is only approximately zero in floating point. + """ + if np.all(simulated == simulated[0]) or np.all(observed == observed[0]): + return None + simulated_deviations = simulated - np.mean(simulated) + observed_deviations = observed - np.mean(observed) + spread = float(np.sqrt(np.sum(simulated_deviations**2) * np.sum(observed_deviations**2))) + return float(np.sum(simulated_deviations * observed_deviations)) / spread + + +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)}, + ) + + +_ENCODING = "utf-8-sig" +"""The encoding both readers open a time series file with. + +It is UTF-8 with the byte order mark removed when the file carries one. The +mark is invisible and belongs to no cell, but it is written at the start of the +first line: left in place it becomes part of the first header cell, where it +turns the time step of a header-less table into a word and hides the refusal +that :func:`_check_csv_header` owes the reader. +""" + +_LAYOUTS = ( + "a time series file is either the table the model writes, whose first line is the " + "header '0;;...' and whose every further line is one time step, or a PCRaster " + "time series, whose header is a title line, the number of columns, and one line per " + "column name starting with the time step column" +) +"""The two layouts a time series file is read in, named by every refusal of a header.""" + + +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. + + :raises ValueError: If the first record is a row of values instead of a + header, which would otherwise read the first time step as the station + ids and lose it from the series. + """ + with path.open(encoding=_ENCODING, 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]] + _check_csv_header(header[0], path) + 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 _check_csv_header(label: str, path: Path) -> None: + """Refuse a ``;``-separated table whose first record is already a time step. + + The first cell of the header is the label of the step column, ``0`` in the + tables the model writes and any non-numeric word in a table written by + hand. A first cell that is another number is a row of values: the file has + no header, its first time step would be read as the station ids and the + remaining columns would be labelled with that step's values. + + :raises ValueError: If the first cell is a number other than ``0``. + """ + try: + step = float(label) + except ValueError: + return + if step == 0: + return + raise ValueError( + f"The time series file {path} has no header: its first line, {label!r}, is " + f"already a time step, and the station ids would be read from it. Add one: " + f"{_LAYOUTS}." + ) + + +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. + + :raises ValueError: If the header is not there. The columns of a time + series carry the ids of the stations they were sampled at, and nothing + in a header-less file says which station a column belongs to; numbering + the columns instead would label them with a guess. + """ + lines = [line for line in path.read_text(encoding=_ENCODING).splitlines() if line.strip()] + if not lines: + raise ValueError(f"The time series file {path} is empty.") + + columns = int(lines[1].strip()) if len(lines) > 1 and lines[1].strip().isdigit() else 0 + names = [line.strip() for line in lines[2 : 2 + columns]] + if columns < 2 or len(names) != columns: + raise ValueError( + f"The time series file {path} has no header: PCRaster time series files carry " + f"one, and its column names are the ids of the stations the columns were " + f"sampled at. Add it: {_LAYOUTS}." + ) + 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 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, 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. Both layouts name their stations in a header, and a + file without one is refused rather than read with numbered columns: the + columns of a time series are stations, and a numbering would compare the + observations of one gauge with the simulation of another. A byte order mark + at the start of the file is dropped with the encoding, so that it cannot + become part of the first header cell. 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, carries no header, 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=_ENCODING) + 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) + + +@dataclass(frozen=True) +class Alignment: + """The steps and the stations two series are compared on. + + :param stations: The station ids both series carry, sorted. + :type stations: list[str] + + :param steps: The time steps both series carry after the spin-up, sorted. + :type steps: numpy.ndarray + + :param simulated_index: Rows of the simulated series, one per shared step. + :type simulated_index: numpy.ndarray + + :param observed_index: Rows of the observed series, one per shared step. + :type observed_index: numpy.ndarray + """ + + stations: list[str] + steps: np.ndarray + simulated_index: np.ndarray + observed_index: np.ndarray + + +def align_series(simulated: Series, observed: Series, spinup_steps: int = 0) -> Alignment: + """Return the steps and the stations two series have in common. + + 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. This is the + one definition of the compared window, used both by the objective and by the + table of the best candidate, so that the two always describe one comparison. + + 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 shared stations, the shared steps and the rows they are at. + :rtype: Alignment + + :raises ValueError: If the two series share no station id or no time step. + """ + 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()}." + ) + + return Alignment( + stations=shared_ids, + steps=shared_steps, + simulated_index=np.flatnonzero(simulated_kept)[simulated_rows], + observed_index=np.flatnonzero(observed_kept)[observed_rows], + ) + + +def evaluate_series( + simulated: Series, + observed: Series, + spinup_steps: int = 0, + stations: Sequence[str] | None = None, +) -> tuple[float, dict[str, StationMetrics]]: + """Return the mean efficiency of a run and the goodness of fit of each station. + + The two series are aligned by :func:`align_series`, and every station they + share is measured, whether or not it enters the objective. ``stations`` + names the ones the mean efficiency is taken over, which is how a run + calibrates on some of its gauges and reports the others as a validation; + without it every shared station enters the mean. The mean is taken over the + selected stations that have an efficiency, so that a station of zero + variance, reported as ``None``, does not drag the mean. + + :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 + + :param stations: Ids of the stations the mean is taken over. ``None``, the + default, takes it over every shared station. + :type stations: collections.abc.Sequence[str], optional + + :return: The mean efficiency of the selected stations and the statistics of + every shared station. + :rtype: tuple[float, dict[str, StationMetrics]] + + :raises ValueError: If the two series share no station id, share no time + step, if the selection names no shared station, or if no selected + station has an efficiency. + """ + alignment = align_series(simulated, observed, spinup_steps) + metrics = { + station: station_metrics( + simulated.stations[station][alignment.simulated_index], + observed.stations[station][alignment.observed_index], + ) + for station in alignment.stations + } + + if stations is None: + selected = list(alignment.stations) + else: + selected = [station for station in alignment.stations if station in set(stations)] + if not selected: + raise ValueError( + f"The station(s) {', '.join(stations) or 'none'} of the objective are " + "none of the stations the two series share " + f"({', '.join(alignment.stations)})." + ) + + values = [metrics[station].nse for station in selected if metrics[station].nse is not None] + if not values: + raise ValueError( + "No station has a Nash-Sutcliffe efficiency: every station of the objective " + f"({', '.join(selected)}) has fewer than two valid pairs or observations " + "without variance." + ) + return float(np.mean(values)), metrics 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..92ebc5f2 --- /dev/null +++ b/tests/unit/calibration/test_objective.py @@ -0,0 +1,695 @@ +"""The calibration objective: masks, Nash-Sutcliffe efficiency, tables and alignment.""" + +import dataclasses + +import numpy as np +import pytest + +from rubem.calibration.objective import ( + INADMISSIBLE_OBJECTIVE, + MISSING_VALUES, + PCRASTER_MISSING, + Series, + StationMetrics, + evaluate_series, + nash_sutcliffe, + objective, + read_series, + station_metrics, + 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, -1.0, 1e31, FLOAT32_MISSING) + + +def nse_of(metrics): + """The efficiency of every station of a mapping of :class:`StationMetrics`.""" + return {station: measured.nse for station, measured in metrics.items()} + + +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("negative", [-1.0, -999.0, -0.5, -9999.0]) + def test_a_negative_value_is_a_gap(self, negative): + # Every quantity the two series carry is a flux or a storage: a + # negative entry is a marker of a gap, whichever one a gauge writes. + assert valid_mask(np.array([negative])).tolist() == [False] + + @pytest.mark.unit + def test_zero_is_a_value_and_not_a_gap(self): + # A station may record no flow at all; only the sign below zero is a gap. + assert valid_mask(np.array([0.0, -0.0])).tolist() == [True, True] + + @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_header_whose_count_carries_trailing_spaces(self, tmp_path): + # The number of columns is read from a line of its own, which a writer + # may pad; the title above it carries spaces of its own and is not read. + path = tmp_path / "observed.tss" + path.write_text( + "timeseries scalar\n3 \ntimestep\n11\n22\n1 1.5 2.5\n2 3.5 4.5\n", + encoding="utf8", + ) + + table = read_series(path) + + assert sorted(table.stations) == ["11", "22"] + assert table.stations["22"].tolist() == [2.5, 4.5] + + @pytest.mark.unit + def test_the_windows_line_endings_of_a_pcraster_time_series(self, tmp_path): + # The carriage return must not be taken for part of the number of + # columns, of a station id or of the last value of a row. + path = tmp_path / "observed.tss" + path.write_bytes( + b"timeseries scalar\r\n3\r\ntimestep\r\n11\r\n22\r\n1 1.5 2.5\r\n2 3.5 4.5\r\n" + ) + + table = read_series(path) + + assert sorted(table.stations) == ["11", "22"] + assert table.steps.tolist() == [1, 2] + assert table.stations["22"].tolist() == [2.5, 4.5] + + @pytest.mark.unit + def test_a_pcraster_time_series_without_the_header_is_refused(self, tmp_path): + # Numbering the columns 1..N would label them with a guess: the columns + # of a time series are stations, and nothing here says which. + path = tmp_path / "observed.tss" + path.write_text("1 1.5 2.5\n2 3.5 4.5\n", encoding="utf8") + + with pytest.raises(ValueError, match="PCRaster time series files carry") as failure: + read_series(path) + + message = str(failure.value) + assert "has no header" in message + assert "0;;..." in message, "the message names both accepted layouts" + assert "title line" in message + + @pytest.mark.unit + @pytest.mark.parametrize( + "content", + [ + "Observed streamflow\n1\ntimestep\n1 1.5\n", + "Observed streamflow\n3\ntimestep\n11\n", + "1 1.5 2.5\n", + ], + ids=["one column", "fewer names than the count", "a single row"], + ) + def test_a_truncated_pcraster_header_is_refused(self, tmp_path, content): + path = tmp_path / "observed.tss" + path.write_text(content, encoding="utf8") + + with pytest.raises(ValueError, match="PCRaster time series files carry"): + read_series(path) + + @pytest.mark.unit + def test_a_csv_table_without_the_header_is_refused(self, tmp_path): + # The first record is already a time step, so the ids would be read + # from it and the first step would be lost from the series. + path = tmp_path / "observed.csv" + path.write_text("1;16.33;22.94\n2;14.86;19.95\n", encoding="utf8") + + with pytest.raises(ValueError, match="has no header") as failure: + read_series(path) + + message = str(failure.value) + assert "'1'" in message, "the message quotes the line it refused" + assert "0;;..." in message, "the message names both accepted layouts" + assert "PCRaster time series" in message + + @pytest.mark.unit + def test_a_byte_order_mark_is_not_part_of_the_step_column_label(self, tmp_path): + # An editor that writes UTF-8 with a byte order mark puts it at the + # start of the first line, where it would otherwise become part of the + # first header cell. + path = tmp_path / "observed.csv" + path.write_bytes("0;A;B\n1;1.5;2.5\n2;3.5;4.5\n".encode("utf8")) + + table = read_series(path) + + assert sorted(table.stations) == ["A", "B"] + assert table.steps.tolist() == [1, 2] + assert table.stations["A"].tolist() == [1.5, 3.5] + + @pytest.mark.unit + def test_a_byte_order_mark_does_not_hide_a_missing_csv_header(self, tmp_path): + # The mark would turn the first time step into a word, and a word is + # the label of a step column: the refusal must not depend on it. + path = tmp_path / "observed.csv" + path.write_bytes("1;16.33;22.94\n2;14.86;19.95\n".encode("utf8")) + + with pytest.raises(ValueError, match="has no header") as failure: + read_series(path) + + assert "'1'" in str(failure.value), "the mark is not part of the line it quotes" + + @pytest.mark.unit + def test_a_byte_order_mark_before_the_title_of_a_time_series(self, tmp_path): + path = tmp_path / "observed.tss" + path.write_bytes( + "timeseries scalar\n3\ntimestep\n11\n22\n1 1.5 2.5\n2 3.5 4.5\n".encode("utf8") + ) + + table = read_series(path) + + assert sorted(table.stations) == ["11", "22"] + assert table.steps.tolist() == [1, 2] + + @pytest.mark.unit + @pytest.mark.parametrize("label", ["0", "timestep", "step", ""]) + def test_the_accepted_spellings_of_the_step_column_label(self, tmp_path, label): + path = tmp_path / "observed.csv" + path.write_text(f"{label};A\n1;1.0\n2;2.0\n", encoding="utf8") + + table = read_series(path) + + assert table.steps.tolist() == [1, 2] + assert table.stations["A"].tolist() == [1.0, 2.0] + + @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 nse_of(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 nse_of(per_station) == {"A": 1.0} + assert per_station["A"].pairs == 2, "only the shared steps are a pair" + 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 nse_of(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 nse_of(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 nse_of(per_station) == {"A": 1.0, "B": None} + assert per_station["B"].pairs == 3, "the station is measured even without an efficiency" + 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) + + +class TestStationMetrics: + @pytest.mark.unit + def test_every_statistic_of_a_hand_computed_pair_of_series(self): + # observed [1, 2, 6], simulated [2, 2, 5]. + # mean_observed = (1 + 2 + 6) / 3 = 3 + # mean_simulated = (2 + 2 + 5) / 3 = 3 + # deviations: observed [-2, -1, 3], simulated [-1, -1, 2] + # std_observed = sqrt((4 + 1 + 9) / 2) = sqrt(7) + # std_simulated = sqrt((1 + 1 + 4) / 2) = sqrt(3) + # r = (2 + 1 + 6) / sqrt(14 * 6) = 9 / sqrt(84) + # rmse = sqrt(((2-1)^2 + 0 + (5-6)^2) / 3) = sqrt(2/3) + # nse = 1 - (1 + 0 + 1) / 14 = 1 - 2/14 = 6/7 + measured = station_metrics(np.array([2.0, 2.0, 5.0]), np.array([1.0, 2.0, 6.0])) + + assert measured.pairs == 3 + assert measured.mean_observed == pytest.approx(3.0) + assert measured.mean_simulated == pytest.approx(3.0) + assert measured.std_observed == pytest.approx(np.sqrt(7.0)) + assert measured.std_simulated == pytest.approx(np.sqrt(3.0)) + assert measured.r == pytest.approx(9.0 / np.sqrt(84.0)) + assert measured.rmse == pytest.approx(np.sqrt(2.0 / 3.0)) + assert measured.nse == pytest.approx(6.0 / 7.0) + + @pytest.mark.unit + def test_it_is_a_frozen_record_of_eight_fields(self): + measured = station_metrics(np.array([1.0, 2.0]), np.array([1.0, 2.0])) + + assert isinstance(measured, StationMetrics) + assert list(dataclasses.asdict(measured)) == [ + "pairs", + "mean_observed", + "std_observed", + "mean_simulated", + "std_simulated", + "r", + "rmse", + "nse", + ] + with pytest.raises(dataclasses.FrozenInstanceError): + measured.nse = 0.0 + + @pytest.mark.unit + def test_only_the_pairs_valid_in_both_series_are_measured(self): + # The third step is a gap in the simulation and the fourth a negative + # observation: the sample is observed [1, 3] and simulated [1, 3]. + # mean 2 on both sides, std sqrt(((1)^2 + (1)^2) / 1) = sqrt(2) + measured = station_metrics( + np.array([1.0, 3.0, np.nan, 9.0]), np.array([1.0, 3.0, 8.0, -1.0]) + ) + + assert measured.pairs == 2 + assert measured.mean_observed == pytest.approx(2.0) + assert measured.std_observed == pytest.approx(np.sqrt(2.0)) + assert measured.rmse == 0.0 + assert measured.nse == 1.0 + + @pytest.mark.unit + def test_a_single_pair_has_a_mean_and_an_error_but_no_spread(self): + # One pair: the mean and the root mean squared error are defined, a + # standard deviation with one degree of freedom and an efficiency are not. + measured = station_metrics(np.array([3.0, np.nan]), np.array([1.0, 2.0])) + + assert measured.pairs == 1 + assert measured.mean_observed == 1.0 + assert measured.mean_simulated == 3.0 + assert measured.rmse == 2.0 + assert measured.std_observed is None + assert measured.std_simulated is None + assert measured.r is None + assert measured.nse is None + + @pytest.mark.unit + def test_a_station_without_one_valid_pair_has_nothing(self): + measured = station_metrics(np.array([1.0, -9999.0]), np.array([np.nan, 2.0])) + + assert measured.pairs == 0 + assert dataclasses.astuple(measured) == (0, None, None, None, None, None, None, None) + + @pytest.mark.unit + def test_a_constant_series_has_no_correlation(self): + # The simulation never moves, so its own spread is zero and the + # correlation has no denominator; the efficiency is still defined, + # since the observations vary. + measured = station_metrics(np.array([2.0, 2.0, 2.0]), np.array([1.0, 2.0, 3.0])) + + assert measured.std_simulated == 0.0 + assert measured.r is None + assert measured.nse == 0.0 + + @pytest.mark.unit + def test_constant_observations_leave_neither_a_correlation_nor_an_efficiency(self): + measured = station_metrics(np.array([0.1, 0.2, 0.3]), np.array([0.1, 0.1, 0.1])) + + # The spread of a constant series is only approximately zero in + # floating point, which is why the guards test the values themselves. + assert measured.std_observed == pytest.approx(0.0, abs=1e-15) + assert measured.r is None + assert measured.nse is None + + @pytest.mark.unit + def test_series_of_different_lengths_are_refused(self): + with pytest.raises(ValueError, match="same shape"): + station_metrics(np.array([1.0, 2.0]), np.array([1.0, 2.0, 3.0])) + + +class TestStationSelection: + @pytest.mark.unit + def test_the_mean_is_taken_over_the_selection_and_every_station_is_measured(self): + # Station A is exact -> 1.0, station B reversed -> -3.0, station C + # exact -> 1.0. The objective averages A and C: (1.0 + 1.0) / 2 = 1.0, + # and B is measured all the same. + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[3.0, 2.0, 1.0], C=[4.0, 5.0, 6.0]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[1.0, 2.0, 3.0], C=[4.0, 5.0, 6.0]) + + mean, per_station = evaluate_series(simulated, observed, stations=("A", "C")) + + assert mean == 1.0 + assert nse_of(per_station) == {"A": 1.0, "B": -3.0, "C": 1.0} + + @pytest.mark.unit + def test_a_selected_station_without_an_efficiency_is_left_out_of_the_mean(self): + # B has constant observations and no efficiency, so the mean of the + # selection is that of A alone. + 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, _ = evaluate_series(simulated, observed, stations=("A", "B")) + + assert mean == 1.0 + + @pytest.mark.unit + def test_a_station_of_the_selection_the_series_do_not_share_is_ignored(self): + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[3.0, 2.0, 1.0]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[1.0, 2.0, 3.0]) + + mean, per_station = evaluate_series(simulated, observed, stations=("A", "Z")) + + assert mean == 1.0 + assert sorted(per_station) == ["A", "B"] + + @pytest.mark.unit + def test_a_selection_that_names_no_shared_station_names_the_ids(self): + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0]) + + with pytest.raises(ValueError, match="none of the stations") as failure: + evaluate_series(simulated, observed, stations=("Y", "Z")) + + message = str(failure.value) + assert "Y, Z" in message + assert "A" in message + + @pytest.mark.unit + def test_no_selected_station_with_an_efficiency_names_the_selection(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, 2.0, 3.0]) + + with pytest.raises(ValueError, match="No station has a Nash-Sutcliffe") as failure: + evaluate_series(simulated, observed, stations=("A",)) + + assert "(A)" in str(failure.value), "only the selection is named" From ffab73e6f411d8e1c774abd3d0b89028d6c62011 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Fri, 18 Sep 2026 16:12:41 -0300 Subject: [PATCH 3/9] Add the decision space of the calibration The search runs over the eight free calibration parameters, with w_3 derived as 1 - w_1 - w_2 and the linear constraint w_1 + w_2 <= 1. decision_space(fixed, bounds) builds the vector of one run: a fixed parameter leaves the vector and is re-inserted in the nine parameters every candidate runs with; a per-run bound may only narrow the range of the application settings; fixing one weight narrows the other so the derived w_3 stays in range and drops the pair constraint, a narrowing that leaves it one value fixes it as well, and fixing both weights checks w_3 before the search. --- rubem/calibration/parameters.py | 511 ++++++++++++++++++++++ tests/unit/calibration/test_parameters.py | 455 +++++++++++++++++++ 2 files changed, 966 insertions(+) create mode 100644 rubem/calibration/parameters.py create mode 100644 tests/unit/calibration/test_parameters.py diff --git a/rubem/calibration/parameters.py b/rubem/calibration/parameters.py new file mode 100644 index 00000000..78c16efd --- /dev/null +++ b/rubem/calibration/parameters.py @@ -0,0 +1,511 @@ +"""The decision vector of the calibration and its translation to model parameters. + +The differential evolution searches the free calibration 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. A +search over all eight is therefore eight-dimensional with the linear constraint +``w_1 + w_2 <= 1``, which keeps the derived weight inside its own range. + +:class:`DecisionSpace`, built by :func:`decision_space`, is that vector for one +run. It narrows the bounds of a parameter to a range of the user's choosing and +takes a parameter out of the vector altogether, at a value the user pins it to: +a fixed parameter keeps its value in the nine parameters every candidate is run +with, but the optimizer never moves it, and the dimension of the search falls by +one. Fixing one of the two searched weights replaces the linear constraint with +a narrower bound on the other, since only one of them is left to vary. + +The bounds a run starts from 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. An override may only narrow them, never widen +them, for the same reason. + +The module-level :func:`bounds`, :func:`vector_to_parameters`, +:func:`parameters_to_vector`, :func:`is_admissible` and +:func:`weights_constraint` are the default space, the one without fixed +parameters and without overridden bounds. + +Nothing here imports SciPy at module level; +:meth:`DecisionSpace.weights_constraint`, the only method that needs it, imports +it when it is called. +""" + +import logging +import math +from collections.abc import Mapping, Sequence +from dataclasses import dataclass +from typing import TYPE_CHECKING, cast + +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.""" + + +@dataclass(frozen=True) +class DecisionSpace: + """The decision vector of one calibration: what is searched and between which bounds. + + The instances are plain data and are built by :func:`decision_space`, which + is what validates them; they travel to the worker processes inside the + evaluation context, so they hold nothing but names and numbers. + + :param free_names: The searched parameters, in the order of the decision + vector: the entries of :data:`FREE_PARAMETERS` that are not fixed. + :type free_names: tuple[str, ...] + + :param fixed: The parameters that are not searched, by name, with the value + every candidate carries. + :type fixed: dict[str, float] + + :param bounds: One ``(minimum, maximum)`` pair per free parameter, in the + order of :attr:`free_names`. + :type bounds: tuple[tuple[float, float], ...] + """ + + free_names: tuple[str, ...] + fixed: dict[str, float] + bounds: tuple[tuple[float, float], ...] + + @property + def dimension(self) -> int: + """Number of coordinates of the decision vector. + + :return: The number of free parameters. + :rtype: int + """ + return len(self.free_names) + + def to_parameters(self, vector: Sequence[float] | np.ndarray) -> dict[str, float]: + """Return the nine calibration parameters a decision vector stands for. + + The free values are read in the order of :attr:`free_names`, the fixed + ones are inserted at the value they were pinned to, 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 free values, in the order of :attr:`free_names`. + :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 :attr:`dimension` entries. + """ + values = np.asarray(vector, dtype=np.float64).ravel() + if values.size != self.dimension: + raise ValueError( + f"The decision vector must have {self.dimension} entries " + f"({', '.join(self.free_names)}), got {values.size}." + ) + parameters = dict(self.fixed) + parameters.update( + {name: float(value) for name, value in zip(self.free_names, 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 from_parameters(self, 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`` and the fixed parameters are ignored, + since they are 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 free values, in the order of :attr:`free_names`. + :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 self.free_names if name not in normalised] + if missing: + raise KeyError( + f"The calibration parameters are missing {', '.join(missing)}; " + f"the decision vector needs {', '.join(self.free_names)}." + ) + return np.asarray([float(normalised[name]) for name in self.free_names], dtype=np.float64) + + def is_admissible(self, 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 weights add up to at most 1 and when the derived ``w_3`` lies + inside its own bound. The derived weight is taken from + :meth:`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 free values, in the order of :attr:`free_names`. + :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 != self.dimension: + return False + if not np.all(np.isfinite(values)): + return False + for value, (minimum, maximum) in zip(values, self.bounds, strict=True): + if not minimum <= value <= maximum: + return False + + parameters = self.to_parameters(values) + if parameters["w_1"] + parameters["w_2"] > 1.0: + return False + minimum, maximum = variable_range(DERIVED_PARAMETER) + return minimum <= parameters[DERIVED_PARAMETER] <= maximum + + def weights_constraint(self) -> "LinearConstraint | None": + """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. The constraint indexes the vector by position, so its + coefficients are placed at the positions the two weights occupy in + :attr:`free_names`. + + A search in which at most one of the two weights is free needs no + constraint: one fixed weight turns the sum into a bound on the other, + which :func:`decision_space` has already applied, and two fixed weights + leave nothing to constrain. + + :return: The constraint over the decision vector, or ``None`` when fewer + than two of the weights are searched. + :rtype: scipy.optimize.LinearConstraint | None + + :raises ImportError: If SciPy is not installed. + """ + positions = [ + self.free_names.index(name) for name in ("w_1", "w_2") if name in self.free_names + ] + if len(positions) < 2: + return None + + from scipy.optimize import LinearConstraint + + coefficients = np.zeros((1, self.dimension), dtype=np.float64) + for position in positions: + coefficients[0, position] = 1.0 + return LinearConstraint(coefficients, -np.inf, 1.0) + + +def decision_space( + fixed: Mapping[str, float] | None = None, + bounds: Mapping[str, tuple[float, float]] | None = None, +) -> DecisionSpace: + """Return the decision space of one calibration. + + Without arguments the space is the whole search: the eight free parameters + of :data:`FREE_PARAMETERS`, each between the bounds of the application + settings. + + ``fixed`` takes parameters out of the decision vector; their values stay in + the nine parameters every candidate is run with. ``bounds`` narrows the range + a free parameter is searched in. Both accept the legacy spellings ``b``, + ``w1`` and ``w2``. + + Fixing one of the two searched weights narrows the upper bound of the other + instead of constraining the pair, so that the derived ``w_3 = 1 - w_1 - w_2`` + cannot fall below its own minimum; when that narrowing leaves the other + weight a single value, it is fixed at that value as well and leaves the + decision vector. Fixing both of them derives ``w_3`` once and checks it here, + before the search starts. + + :param fixed: The parameters that are not searched, with the value every + candidate carries. Defaults to ``None``, no fixed parameter. + :type fixed: collections.abc.Mapping[str, float], optional + + :param bounds: The ``(minimum, maximum)`` range to search a parameter in, + by name. Defaults to ``None``, the ranges of the application settings. + :type bounds: collections.abc.Mapping[str, tuple[float, float]], optional + + :return: The decision vector of the run. + :rtype: DecisionSpace + + :raises ValueError: If a name is not a searched parameter or is the derived + ``w_3``, if a fixed value lies outside the range of the settings, if an + overridden bound is not a narrower range of the settings range, if a + fixed weight leaves the other one no admissible value, if two fixed + weights derive a ``w_3`` outside its range, or if every parameter is + fixed and nothing is left to search. + """ + fixed_values = { + _canonical(name, "fixed"): float(value) for name, value in (fixed or {}).items() + } + overrides = {_canonical(name, "bounded"): value for name, value in (bounds or {}).items()} + + settings_ranges = {name: variable_range(name) for name in FREE_PARAMETERS} + effective = dict(settings_ranges) + for name, override in overrides.items(): + if name in fixed_values: + raise ValueError( + f"'{name}' is fixed at {fixed_values[name]}, so a bound for it has nowhere " + "to apply: only the free parameters are searched." + ) + effective[name] = _narrowed(name, override, settings_ranges[name]) + + for name, value in fixed_values.items(): + minimum, maximum = settings_ranges[name] + if not math.isfinite(value) or not minimum <= value <= maximum: + raise ValueError( + f"'{name}' cannot be fixed at {value}: the value lies outside the range " + f"({minimum}, {maximum}) of the application settings." + ) + + free_names = tuple(name for name in FREE_PARAMETERS if name not in fixed_values) + if not free_names: + raise ValueError( + "Every calibration parameter is fixed, so the search has nothing to look for; " + "leave at least one of them free." + ) + + effective.update(_weight_bounds(fixed_values, effective)) + # A fixed weight can leave the other one a single admissible value (w_1 = 1 + # leaves w_2 = 0). A parameter with one value is not searched: it joins the + # fixed ones, so the dimension, the budget and the initial population describe + # the search that really runs, and no coordinate of zero width reaches SciPy. + for name in ("w_1", "w_2"): + if name in fixed_values: + continue + minimum, maximum = effective[name] + if minimum == maximum: + fixed_values[name] = minimum + free_names = tuple(free for free in free_names if free != name) + if not free_names: + raise ValueError( + "Every calibration parameter is fixed, so the search has nothing to look for; " + "leave at least one of them free." + ) + return DecisionSpace( + free_names=free_names, + fixed=fixed_values, + bounds=tuple(effective[name] for name in free_names), + ) + + +def _canonical(name: str, action: str) -> str: + """Return the canonical spelling of a parameter the caller named. + + :raises ValueError: If the name is the derived weight or is not a searched + parameter at all. + """ + canonical = _ALIASES.get(name, name) + if canonical == DERIVED_PARAMETER: + raise ValueError( + f"'{name}' cannot be {action}: the slope factor weight is derived from the " + "other two weights (w_3 = 1 - w_1 - w_2) and is not part of the decision vector; " + "fix or bound w_1 and w_2 instead." + ) + if canonical not in FREE_PARAMETERS: + raise ValueError( + f"'{name}' cannot be {action}: it is not a searched calibration parameter. " + f"The searched parameters are {', '.join(FREE_PARAMETERS)}." + ) + return canonical + + +def _narrowed( + name: str, override: Sequence[float], settings_range: tuple[float, float] +) -> tuple[float, float]: + """Return an overridden bound, checked against the range of the settings. + + :raises ValueError: If the override is not a range with a minimum below its + maximum, or if it is not contained in the range of the settings. + """ + try: + values = tuple(float(value) for value in override) + except TypeError: + # A single number instead of a pair: the same mistake as a pair of the + # wrong length, and it is answered with the same sentence. + raise ValueError( + f"The bound of '{name}' must be a (minimum, maximum) pair, got {override!r}." + ) from None + if len(values) != 2: + raise ValueError( + f"The bound of '{name}' must be a (minimum, maximum) pair, got {len(values)} value(s)." + ) + minimum, maximum = values + settings_minimum, settings_maximum = settings_range + if ( + not math.isfinite(minimum) + or not math.isfinite(maximum) + or minimum >= maximum + or minimum < settings_minimum + or maximum > settings_maximum + ): + raise ValueError( + f"The bound ({minimum}, {maximum}) of '{name}' is not a narrower range of its " + f"application settings range ({settings_minimum}, {settings_maximum}): a bound " + "narrows the settings range and its minimum lies below its maximum." + ) + return minimum, maximum + + +def _weight_bounds( + fixed_values: Mapping[str, float], effective: Mapping[str, tuple[float, float]] +) -> dict[str, tuple[float, float]]: + """Return the bounds the fixed weights impose on the weight that is left. + + With both weights searched the pair is handled by the linear constraint of + :meth:`DecisionSpace.weights_constraint` and nothing is narrowed here. With + one of them fixed the constraint would have a single term, which is a bound: + the free weight may not exceed what leaves the derived ``w_3`` at or above + its own minimum. With both of them fixed the derived weight is a number, and + it is checked here. + + :raises ValueError: If the free weight has no admissible value left, or if + two fixed weights derive a ``w_3`` outside its range. + """ + derived_minimum, derived_maximum = variable_range(DERIVED_PARAMETER) + weights = [name for name in ("w_1", "w_2") if name in fixed_values] + if not weights: + return {} + + if len(weights) == 2: + derived = 1.0 - (fixed_values["w_1"] + fixed_values["w_2"]) + if not derived_minimum <= derived <= derived_maximum: + raise ValueError( + f"Fixing w_1 at {fixed_values['w_1']} and w_2 at {fixed_values['w_2']} derives " + f"w_3 = {derived}, outside its range ({derived_minimum}, {derived_maximum}); " + "the three weights must add up to 1." + ) + return {} + + (pinned,) = weights + free_weight = "w_2" if pinned == "w_1" else "w_1" + minimum, maximum = effective[free_weight] + narrowed = min(maximum, 1.0 - fixed_values[pinned] - derived_minimum) + if narrowed < minimum: + raise ValueError( + f"Fixing {pinned} at {fixed_values[pinned]} leaves {free_weight} no value in its " + f"range ({minimum}, {maximum}): the derived w_3 = 1 - w_1 - w_2 would fall below " + f"{derived_minimum} for every one of them." + ) + return {free_weight: (minimum, narrowed)} + + +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 of a + search without fixed parameters and without overridden bounds. + + :return: One ``(minimum, maximum)`` pair per free parameter. + :rtype: list[tuple[float, float]] + """ + return list(decision_space().bounds) + + +def vector_to_parameters(vector: Sequence[float] | np.ndarray) -> dict[str, float]: + """Return the nine calibration parameters a decision vector stands for. + + The vector is the one of the default decision space, the eight free values + in the order of :data:`FREE_PARAMETERS`; see + :meth:`DecisionSpace.to_parameters`. + + :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. + """ + return decision_space().to_parameters(vector) + + +def parameters_to_vector(parameters: Mapping[str, float]) -> np.ndarray: + """Return the decision vector of a set of calibration parameters. + + The vector is the one of the default decision space; see + :meth:`DecisionSpace.from_parameters`. + + :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. + """ + return decision_space().from_parameters(parameters) + + +def is_admissible(vector: Sequence[float] | np.ndarray) -> bool: + """Whether a decision vector describes a configuration the model accepts. + + The vector is the one of the default decision space; see + :meth:`DecisionSpace.is_admissible`. + + :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 + """ + return decision_space().is_admissible(vector) + + +def weights_constraint() -> "LinearConstraint": + """Return the linear constraint ``w_1 + w_2 <= 1`` on the decision vector. + + The constraint is the one of the default decision space, in which both + weights are searched; see :meth:`DecisionSpace.weights_constraint`. + + :return: The constraint, over the eight entries of the decision vector. + :rtype: scipy.optimize.LinearConstraint + + :raises ImportError: If SciPy is not installed. + """ + # Both weights are searched in the default space, so the constraint is there. + return cast("LinearConstraint", decision_space().weights_constraint()) diff --git a/tests/unit/calibration/test_parameters.py b/tests/unit/calibration/test_parameters.py new file mode 100644 index 00000000..a2c209a7 --- /dev/null +++ b/tests/unit/calibration/test_parameters.py @@ -0,0 +1,455 @@ +"""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, + decision_space, + 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]) + + +class TestDefaultDecisionSpace: + @pytest.mark.unit + def test_a_space_without_arguments_is_the_whole_search(self): + space = decision_space() + + assert space.free_names == FREE_PARAMETERS + assert space.fixed == {} + assert space.dimension == 8 + assert list(space.bounds) == bounds() + + @pytest.mark.unit + def test_the_module_level_functions_are_the_default_space(self): + space = decision_space() + vector = admissible_vector() + parameters = vector_to_parameters(vector) + + assert space.to_parameters(vector) == parameters + assert space.from_parameters(parameters) == pytest.approx(parameters_to_vector(parameters)) + assert space.is_admissible(vector) is is_admissible(vector) is True + assert np.array_equal( + np.asarray(space.weights_constraint().A), np.asarray(weights_constraint().A) + ) + + +class TestFixedParameters: + @pytest.mark.unit + def test_a_fixed_parameter_leaves_the_vector_and_stays_in_the_parameters(self): + space = decision_space(fixed={"x": 0.0}) + + assert space.free_names == tuple(name for name in FREE_PARAMETERS if name != "x") + assert space.dimension == 7 + assert space.fixed == {"x": 0.0} + assert list(space.bounds) == [ + SETTINGS_BOUNDS[name] for name in FREE_PARAMETERS if name != "x" + ] + + parameters = space.to_parameters(np.delete(admissible_vector(), FREE_PARAMETERS.index("x"))) + + assert list(parameters) == list(CALIBRATION_PARAMETERS) + assert parameters["x"] == 0.0 + assert parameters["alpha"] == 4.5 + assert parameters["w_3"] == pytest.approx(1.0 - 0.666) + + @pytest.mark.unit + def test_the_free_vector_of_a_space_with_fixed_parameters_round_trips(self): + space = decision_space(fixed={"x": 0.25, "alpha": 2.0}) + parameters = vector_to_parameters(admissible_vector()) + + vector = space.from_parameters(parameters) + + assert len(vector) == space.dimension == 6 + # The fixed coordinates are dropped on the way in and restored, at the + # value they were pinned to, on the way out. + assert space.to_parameters(vector) == pytest.approx( + {**parameters, "x": 0.25, "alpha": 2.0}, abs=1e-12 + ) + + @pytest.mark.unit + def test_the_legacy_spellings_name_the_parameter_that_is_fixed(self): + space = decision_space(fixed={"b": 0.25, "w1": 0.5}) + + assert space.fixed == {"beta": 0.25, "w_1": 0.5} + assert "beta" not in space.free_names + + @pytest.mark.unit + def test_a_vector_of_the_wrong_size_names_the_free_parameters(self): + space = decision_space(fixed={"x": 0.0}) + + with pytest.raises(ValueError, match="must have 7 entries"): + space.to_parameters(admissible_vector()) + + @pytest.mark.unit + def test_the_derived_weight_cannot_be_fixed(self): + with pytest.raises(ValueError, match="w_3"): + decision_space(fixed={"w_3": 0.3}) + with pytest.raises(ValueError, match="derived"): + decision_space(fixed={"w3": 0.3}) + + @pytest.mark.unit + def test_an_unknown_parameter_cannot_be_fixed(self): + with pytest.raises(ValueError, match="not a searched calibration parameter"): + decision_space(fixed={"gamma": 1.0}) + + @pytest.mark.unit + @pytest.mark.parametrize("value", [-0.5, 1.5, float("nan")]) + def test_a_fixed_value_outside_the_settings_range_is_refused(self, value): + with pytest.raises(ValueError, match=r"'x' cannot be fixed"): + decision_space(fixed={"x": value}) + + @pytest.mark.unit + def test_fixing_every_parameter_leaves_nothing_to_search(self): + every = {name: minimum for name, (minimum, _) in SETTINGS_BOUNDS.items()} + every.pop(DERIVED_PARAMETER) + + with pytest.raises(ValueError, match="nothing to look for"): + decision_space(fixed=every) + + +class TestFixedWeights: + @pytest.mark.unit + def test_a_weight_left_with_one_value_is_fixed_as_well(self): + """w_1 = 1 leaves w_2 no room: w_2 = 0 is fixed and leaves the vector. + + The dimension, the budget and the initial population then describe the + search that runs, and no zero-width coordinate reaches the optimizer. + """ + space = decision_space(fixed={"w_1": 1.0}) + + assert "w_2" not in space.free_names + assert space.fixed == {"w_1": 1.0, "w_2": 0.0} + assert space.dimension == 6 + assert space.weights_constraint() is None + parameters = space.to_parameters([4.5, 0.5, 5.0, 0.5, 0.5, 0.5]) + assert (parameters["w_1"], parameters["w_2"], parameters["w_3"]) == (1.0, 0.0, 0.0) + + @pytest.mark.unit + def test_fixing_one_weight_narrows_the_other_instead_of_constraining_the_pair(self): + space = decision_space(fixed={"w_1": 0.7}) + position = space.free_names.index("w_2") + + # w_3 = 1 - w_1 - w_2 may not fall below its own minimum, so w_2 is + # searched in [0, 0.3] and the pair needs no constraint of its own. + assert space.bounds[position] == (0.0, pytest.approx(0.3)) + assert space.weights_constraint() is None + + @pytest.mark.unit + def test_the_narrowed_bound_keeps_the_override_when_the_override_is_tighter(self): + space = decision_space(fixed={"w_2": 0.5}, bounds={"w_1": (0.1, 0.2)}) + position = space.free_names.index("w_1") + + assert space.bounds[position] == (0.1, 0.2) + + @pytest.mark.unit + def test_a_fixed_weight_that_leaves_the_other_no_value_is_refused(self): + with pytest.raises(ValueError, match="no value in its range"): + decision_space(fixed={"w_1": 0.9}, bounds={"w_2": (0.2, 0.5)}) + + @pytest.mark.unit + def test_fixing_both_weights_checks_the_derived_one(self): + space = decision_space(fixed={"w_1": 0.3, "w_2": 0.4}) + vector = np.array([4.5, 0.5, 5.0, 0.5, 0.5, 0.5], dtype=np.float64) + + assert space.free_names == ("alpha", "beta", "rcd", "f", "alpha_gw", "x") + assert space.weights_constraint() is None + assert space.to_parameters(vector)["w_3"] == pytest.approx(0.3) + assert space.is_admissible(vector) + + @pytest.mark.unit + def test_two_fixed_weights_that_add_up_above_one_are_refused(self): + with pytest.raises(ValueError, match="derives w_3"): + decision_space(fixed={"w_1": 0.6, "w_2": 0.6}) + + @pytest.mark.unit + def test_the_constraint_follows_the_positions_of_the_weights_in_the_vector(self): + pytest.importorskip("scipy.optimize") + space = decision_space(fixed={"alpha": 5.0}) + constraint = space.weights_constraint() + + expected = np.zeros((1, space.dimension)) + expected[0, space.free_names.index("w_1")] = 1.0 + expected[0, space.free_names.index("w_2")] = 1.0 + assert np.array_equal(np.asarray(constraint.A), expected) + assert np.all(np.asarray(constraint.ub) == 1.0) + + +class TestOverriddenBounds: + @pytest.mark.unit + def test_an_override_narrows_the_range_of_the_settings(self): + space = decision_space(bounds={"rcd": (2.0, 5.0)}) + position = space.free_names.index("rcd") + + assert space.bounds[position] == (2.0, 5.0) + assert [ + bound + for name, bound in zip(space.free_names, space.bounds, strict=True) + if name != "rcd" + ] == [SETTINGS_BOUNDS[name] for name in FREE_PARAMETERS if name != "rcd"] + + @pytest.mark.unit + def test_a_candidate_outside_the_narrowed_bound_is_not_admissible(self): + space = decision_space(bounds={"rcd": (2.0, 5.0)}) + vector = admissible_vector() + + assert space.is_admissible(vector) + vector[FREE_PARAMETERS.index("rcd")] = 8.0 + assert not space.is_admissible(vector) + # The same candidate is admissible in the default space: the refusal is + # the narrowed bound of this run and not a bound of the model. + assert is_admissible(vector) + + @pytest.mark.unit + @pytest.mark.parametrize( + "override", + [(0.5, 5.0), (2.0, 12.0), (5.0, 2.0), (3.0, 3.0), (float("nan"), 5.0)], + ids=["below", "above", "reversed", "empty", "not a number"], + ) + def test_a_bound_that_does_not_narrow_the_settings_range_is_refused(self, override): + with pytest.raises(ValueError, match=r"of 'rcd' is not a narrower range") as failure: + decision_space(bounds={"rcd": override}) + + # Both ranges are named, so the message says what was asked and what is + # available. + assert "(1.0, 10.0)" in str(failure.value) + + @pytest.mark.unit + @pytest.mark.parametrize( + "override", [3.0, (3.0,), (1.0, 2.0, 3.0)], ids=["a number", "one value", "three values"] + ) + def test_a_bound_that_is_not_a_pair_is_refused(self, override): + # A bound is a range: anything else is answered with the pair it needs + # and not with the error of the conversion it failed. + with pytest.raises(ValueError, match=r"bound of 'rcd' must be a \(minimum, maximum\) pair"): + decision_space(bounds={"rcd": override}) + + @pytest.mark.unit + def test_a_bound_of_the_derived_or_of_an_unknown_parameter_is_refused(self): + with pytest.raises(ValueError, match="derived"): + decision_space(bounds={"w_3": (0.1, 0.5)}) + with pytest.raises(ValueError, match="not a searched calibration parameter"): + decision_space(bounds={"gamma": (0.1, 0.5)}) + + @pytest.mark.unit + def test_a_bound_of_a_fixed_parameter_has_nowhere_to_apply(self): + with pytest.raises(ValueError, match="nowhere to apply"): + decision_space(fixed={"rcd": 3.0}, bounds={"rcd": (2.0, 5.0)}) From da951c514a2bb80819edbf9c2c54e6c29ebd6379 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Fri, 18 Sep 2026 16:12:41 -0300 Subject: [PATCH 4/9] Add the differential evolution calibration runner _worker.evaluate is the module-level entry point a spawned worker runs: an inadmissible candidate is rejected without a run; otherwise the format 1.0 document of the configuration is copied with the candidate's parameters, every raster series disabled and only the target time series enabled into a temporary directory, the model runs in-process through the Python API, the table is compared with the observed series and one JSON record (parameters, NSE, per-station metrics, objective, start time, error) is written atomically; failures are recorded and ranked last. runner.calibrate loads and validates the configuration once, refuses an observed series that shares no time step or station with what the configuration samples, writes observed.csv, and drives SciPy's differential_evolution over a spawn ProcessPoolExecutor with one task per silenced worker. It then consolidates evaluations.csv, runs the best candidate once more to write best_.csv and stations.csv, and writes result.json and the calibrated configuration in the input's own format. Progress goes through the rubem.progress logger; an interrupted or failed search still consolidates the evaluations, and a worker killed by the system is reported with the advice to lower the worker count. --- rubem/calibration/_worker.py | 384 +++++++ rubem/calibration/runner.py | 1379 +++++++++++++++++++++++++ tests/unit/calibration/test_runner.py | 1233 ++++++++++++++++++++++ tests/unit/calibration/test_worker.py | 383 +++++++ 4 files changed, 3379 insertions(+) create mode 100644 rubem/calibration/_worker.py create mode 100644 rubem/calibration/runner.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..b0d1e8b9 --- /dev/null +++ b/rubem/calibration/_worker.py @@ -0,0 +1,384 @@ +"""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. + +:func:`evaluate` never 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. :func:`simulate_best`, which the search +calls once at the end to keep the series of the candidate it chose, does raise: +there is one named candidate and nothing to rank a failure behind. +""" + +import copy +import dataclasses +import json +import logging +import os +import shutil +import tempfile +import time +import uuid +from collections.abc import Sequence +from dataclasses import dataclass, field +from datetime import UTC, datetime +from pathlib import Path + +import numpy as np + +from ..configuration.model_configuration_file_v1 import VARIABLE_IDS +from .objective import ( + INADMISSIBLE_OBJECTIVE, + Series, + StationMetrics, + evaluate_series, + objective, + read_series, +) +from .parameters import DecisionSpace, decision_space + +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 + + :param space: The decision vector of the run, which turns a candidate into + the nine parameters of the model. Defaults to the whole search, without + fixed parameters and without narrowed bounds. + :type space: rubem.calibration.parameters.DecisionSpace + + :param stations: Ids of the stations the objective averages. ``None``, the + default, averages every station the two series share; the others are + measured either way. + :type stations: tuple[str, ...] | None + """ + + document: dict + base_dir: str | None + variable: str + observed: Series + spinup_steps: int + temp_dir: str + evaluations_dir: str + validate_input: bool = False + space: DecisionSpace = field(default_factory=decision_space) + stations: tuple[str, ...] | None = None + + +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 free values in the order of the + :attr:`~rubem.calibration.parameters.DecisionSpace.free_names` of the + decision space of the context. + :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_at = datetime.now(UTC).isoformat() + started = time.perf_counter() + try: + parameters = context.space.to_parameters(vector) + except ValueError as error: + _write_record( + context, started_at, {}, None, {}, INADMISSIBLE_OBJECTIVE, 0.0, _describe(error) + ) + return INADMISSIBLE_OBJECTIVE + + if not context.space.is_admissible(vector): + logger.debug("Rejecting an inadmissible candidate without running the model.") + _write_record( + context, + started_at, + 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, metrics = _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, + started_at, + 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, + started_at, + parameters, + nse, + metrics, + value, + time.perf_counter() - started, + None, + ) + return value + + +def simulate_best( + document: dict, + base_dir: str | None, + parameters: dict[str, float], + output_dir: str, + variable: str, +) -> str: + """Run one candidate and keep what it wrote. + + This is how the calibration obtains the series of the candidate it chose: + the run is the very run an evaluation makes, from the same derived + configuration, but the output directory is left in place so that the table + can be read afterwards. It lives at module level and takes plain data for + the same reason :func:`evaluate` does: the spawned workers import this + module by name to find the function they have to call. + + Unlike :func:`evaluate` it raises: the caller asked for one named + candidate's series and has nothing to rank a failure behind. + + :param document: The calibrated configuration as a format 1.0 document. + :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 parameters: The nine calibration parameters of the candidate. + :type parameters: dict[str, float] + + :param output_dir: Directory the run writes into. It is not removed. + :type output_dir: str + + :param variable: Id of the output variable the table is wanted for. + :type variable: str + + :return: The path of the CSV table the run wrote for ``variable``. + :rtype: str + + :raises RuntimeError: If the run wrote no CSV table for the variable. + """ + from ..api import Model + + derived = _document_for(document, variable, parameters, output_dir) + result = Model.from_config(derived, validate_input=False, base_dir=base_dir).run() + return str(_time_series_table(result, variable)) + + +def _run_and_evaluate( + context: EvaluationContext, + parameters: dict[str, float], + output_dir: str, +) -> tuple[float, dict[str, StationMetrics]]: + """Run the candidate and return its mean efficiency and the statistics 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, context.stations) + + +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 of an evaluation.""" + return _document_for(context.document, context.variable, parameters, output_dir) + + +def _document_for( + base_document: dict, + variable: str, + 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 + directory of the run, 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(base_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), + 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, + started_at: str, + parameters: dict[str, float], + nse: float | None, + metrics: dict[str, StationMetrics], + 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. + + The record carries the goodness of fit of every station under + ``station_metrics`` and, under ``station_nse``, the efficiencies alone, + taken from the very same statistics so that the two can never disagree. + ``started_at`` is the wall-clock moment the evaluation began, in UTC: the + records are consolidated in the order of the candidates they evaluated, and + this is what puts them back in the order they were made in. + """ + station_metrics = { + station: dataclasses.asdict(measured) for station, measured in metrics.items() + } + record = { + "id": uuid.uuid4().hex, + "pid": os.getpid(), + "started_at": started_at, + "parameters": parameters, + "nse": nse, + "station_nse": {station: measured["nse"] for station, measured in station_metrics.items()}, + "station_metrics": station_metrics, + "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/runner.py b/rubem/calibration/runner.py new file mode 100644 index 00000000..8ffc5501 --- /dev/null +++ b/rubem/calibration/runner.py @@ -0,0 +1,1379 @@ +"""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, before the search, ``observed.csv``, the coverage +and the summary statistics of the observations over the compared steps, and +after it ``evaluations.csv``, one row per evaluation, ``stations.csv``, the +goodness of fit of the best candidate at every station, ``best_.csv``, +the observed and the simulated series of that candidate side by side, +``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 of evaluations are ordered by the +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; the ``started_at`` column is what puts them back in the order they were +made in. It does not make the table itself reproducible: the ``id``, ``pid``, +``started_at`` 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 * N)``, + with ``N`` the number of free parameters, 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 and the eight free parameters of a search + without fixed ones 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 shutil +import sys +import tempfile +from concurrent.futures import ProcessPoolExecutor +from concurrent.futures.process import BrokenProcessPool +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._ranges import variable_range +from ..configuration.model_configuration_file_v1 import ( + VARIABLE_IDS, + Aggregation, + ModelConfigurationFileV1, +) +from ._worker import _INADMISSIBLE_ERROR, EvaluationContext, evaluate, simulate_best +from .objective import Series, align_series, read_series, valid_mask +from .parameters import ( + CALIBRATION_PARAMETERS, + FREE_PARAMETERS, + DecisionSpace, + decision_space, +) + +logger = logging.getLogger(__name__) + +progress = logging.getLogger("rubem.progress") +"""The logger of what a person watching a calibration wants to see. + +The budget, one line per generation and the closing summary go here; the module +logger keeps the diagnostics, the warnings and the refusals. The command line +routes this logger to standard output, as it does for a simulation, so that an +embedded calibration stays silent unless its host asks for the progress. +""" + +EVALUATIONS_DIRNAME = "evaluations" +"""Name of the directory under the run directory that collects the JSON records.""" + +EVALUATIONS_CSV = "evaluations.csv" +RESULT_JSON = "result.json" +OBSERVED_CSV = "observed.csv" +STATIONS_CSV = "stations.csv" + +EVALUATION_COLUMNS = ( + "id", + "pid", + "started_at", + *CALIBRATION_PARAMETERS, + "nse", + "objective", + "elapsed_seconds", + "error", +) +"""The columns of ``evaluations.csv``, one row per evaluation.""" + +OBSERVED_COLUMNS = ( + "station", + "in_selection", + "pairs_in_window", + "dropped", + "mean", + "std", + "min", + "max", +) +"""The columns of ``observed.csv``, one row per station of the observed series.""" + +STATION_COLUMNS = ( + "station", + "in_selection", + "pairs", + "mean_observed", + "std_observed", + "mean_simulated", + "std_simulated", + "r", + "rmse", + "nse", +) +"""The columns of ``stations.csv``, one row per station the best candidate compared.""" + +STATION_SEPARATOR = ";" +"""The separator of the per-station tables, the one the model writes its own tables with.""" + + +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, N)`` array with ``S > 4`` + members and one column per free parameter, which is one column fewer + for every parameter :attr:`fixed` takes out of the decision vector. + :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 + + :param bounds: The ``(minimum, maximum)`` range to search a parameter in, by + name; it narrows the range of the application settings, which is what a + parameter left out of the mapping is searched in. ``None``, the default, + searches every parameter in the range of the settings. + :type bounds: dict[str, tuple[float, float]] | None + + :param fixed: The parameters that are not searched, by name, with the value + every candidate carries. A fixed parameter leaves the decision vector, + so the search loses one dimension, and keeps its value in the nine + parameters of every run and of the result. ``None``, the default, + searches every parameter. + :type fixed: dict[str, float] | None + + :param stations: Ids of the stations whose efficiency the objective averages. + ``None``, the default, uses every station the observed series and the + run have in common. + :type stations: tuple[str, ...] | None + """ + + 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 + bounds: dict[str, tuple[float, float]] | None = None + fixed: dict[str, float] | None = None + stations: tuple[str, ...] | None = None + + +@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 + + :param best_series: The observed and the simulated series of the best + candidate, side by side on the compared steps. + :type best_series: pathlib.Path + + :param stations_csv: The goodness of fit of the best candidate at every + station the two series share. + :type stations_csv: 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 + best_series: Path + stations_csv: 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 and efficiency so far and the records. + + The best efficiency is read from the records the workers have written, + which is where the efficiencies live: what SciPy reports is the + objective, and the efficiency behind it is what a person watching a + calibration reads. + + 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) + records = _read_records(self.evaluations_dir) + best = _best_record(records) + best_nse = None if best is None else best.get("nse") + progress.info( + "Generation %d: best objective %.6g, best NSE %s, %d evaluation(s) recorded.", + self.generations, + float(intermediate_result.fun), + "n/a" if best_nse is None else f"{float(best_nse):.6f}", + len(records), + ) + 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, if the observed series shares no time + step or no station with what the configuration will sample, if the + station selection names none of those stations, if a worker process + died before finishing its evaluation, or if every evaluation of the run + failed. An interrupted or crashed search still leaves the table of the + evaluations it made behind. + :raises ValueError: If the decision space the settings describe is not a + search (an unknown or derived parameter, a fixed value or an overridden + bound outside the range of the application settings, every parameter + fixed), or if an explicit initial population does not have one column per + free parameter. + :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}." + ) + _check_stations(settings.stations) + + space = decision_space(fixed=settings.fixed, bounds=settings.bounds) + if not isinstance(settings.init, str): + shape = np.shape(settings.init) + if len(shape) != 2 or shape[1] != space.dimension: + raise ValueError( + f"The initial population has the shape {shape} and the search has " + f"{space.dimension} free parameter(s) ({', '.join(space.free_names)}); it must " + f"be a two-dimensional array of {shape[0] if shape else 0} member(s) with one " + "column per free parameter." + ) + + 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) + aggregation = file_v1.model_simulation_output.time_series_samples.aggregation + sampled_ids = ( + None + if aggregation is Aggregation.ZONES + else _station_ids(configuration.raster_files.sample_locations) + ) + compared_steps = _compared_steps( + first_step=configuration.simulation_period.first_step, + last_step=configuration.simulation_period.last_step, + spinup_steps=settings.spinup_steps, + ) + _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=sampled_ids, + ) + _check_station_selection(settings.stations, observed, sampled_ids) + x0 = _starting_point(space, 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." + ) + observed_csv = run_directory / OBSERVED_CSV + _write_observed_csv(observed_csv, observed, compared_steps, settings.stations) + _report_observed_coverage(observed, compared_steps, settings.stations) + 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, + space=space, + stations=settings.stations, + ) + + planned_members = _population_size(settings, space.dimension) + progress.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).", + space.dimension, + planned_members, + settings.maxiter, + planned_members * (settings.maxiter + 1), + workers, + ) + _report_decision_space(space, settings.stations) + + # A search with at most one free weight has no constraint to hand over: the + # sum of the weights is then a bound of the decision space, and SciPy takes + # a different code path for a constrained search, so the argument is left + # out instead of being passed as ``None``. + constraint = space.weights_constraint() + callback = _Progress(evaluations_dir) + executor = _pool(workers) + arguments: dict[str, Any] = { + "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, + "callback": callback, + } + if constraint is not None: + arguments["constraints"] = constraint + evaluations_csv = run_directory / EVALUATIONS_CSV + try: + try: + optimum = differential_evolution( + functools.partial(evaluate, context=context), + list(space.bounds), + **arguments, + ) + finally: + executor.shutdown(wait=True) + except BaseException as error: + # An interruption or a crash after the first evaluation still leaves a + # table: the evaluations that were made are hours of model runs, and a + # calibration that has to be stopped is read from them like any other. + _consolidate(evaluations_csv, evaluations_dir) + if isinstance(error, BrokenProcessPool): + raise CalibrationError( + "A worker process died before finishing its evaluation, most often " + "because the system killed it for lack of memory: every worker holds " + "the rasters of one model run. Run the calibration again with a lower " + f"--workers. The evaluations made so far are in {evaluations_csv}." + ) from error + raise + + records = _read_records(evaluations_dir) + _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 = space.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" + stations_csv = run_directory / STATIONS_CSV + best_series = run_directory / f"best_{settings.variable}.csv" + _write_stations_csv(stations_csv, best.get("station_metrics") or {}, settings.stations) + _write_calibrated_config(calibrated_config, configuration, file_v1, best_parameters) + _write_best_series( + best_series, + document=file_v1.to_dict(), + base_dir=configuration.base_dir, + parameters=best_parameters, + variable=settings.variable, + observed=observed, + spinup_steps=settings.spinup_steps, + temp_dir=temp_dir, + workers=workers, + evaluations_csv=evaluations_csv, + ) + _write_result_json( + result_json, + best_parameters=best_parameters, + best_nse=best["nse"], + best_objective=best_objective, + optimum=optimum, + settings=settings, + space=space, + workers=workers, + temp_dir=temp_dir, + population_size=population_size, + artifacts={ + "evaluations_csv": str(evaluations_csv), + "observed_csv": str(observed_csv), + "stations_csv": str(stations_csv), + "best_series": str(best_series), + "calibrated_config": str(calibrated_config), + }, + ) + + progress.info( + "Calibration finished after %d evaluation(s) in %d generation(s): best objective " + "%.6g, best NSE %s.", + int(optimum.nfev), + int(optimum.nit), + best_objective, + "n/a" if best["nse"] is None else f"{float(best['nse']):.6f}", + ) + 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, + best_series=best_series, + stations_csv=stations_csv, + ) + + +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 _check_stations(stations: tuple[str, ...] | None) -> None: + """Refuse a station selection the objective could not be averaged over. + + :param stations: The ids the caller asked the objective to average, or + ``None`` for every station the two series share. + :type stations: tuple[str, ...] | None + + :raises CalibrationError: If the selection is not a non-empty sequence of + non-empty station ids. + """ + if stations is None: + return + if isinstance(stations, str) or not isinstance(stations, (tuple, list)): + raise CalibrationError( + "The stations of the objective must be given as a tuple of station ids, " + f"got {type(stations).__name__}." + ) + if not stations: + raise CalibrationError( + "The stations of the objective are an empty selection; name at least one station, " + "or leave them unset to average every station the two series share." + ) + wrong = [station for station in stations if not isinstance(station, str) or not station.strip()] + if wrong: + raise CalibrationError( + f"The station id(s) {', '.join(repr(station) for station in wrong)} are not " + "non-empty station ids; a station is named by the id it carries in the header of " + "the time series table." + ) + + +def _pool(workers: int) -> ProcessPoolExecutor: + """Return the process pool one evaluation at a time is run in. + + ``spawn`` with one task per worker is what isolates the runs from each + other: PCRaster keeps its clone and its raster memory process-wide, so a + worker reused for a second evaluation would carry the state of the first. + """ + return ProcessPoolExecutor( + max_workers=workers, + mp_context=multiprocessing.get_context("spawn"), + max_tasks_per_child=1, + initializer=_prepare_worker, + ) + + +def _consolidate(evaluations_csv: Path, evaluations_dir: Path) -> None: + """Write the table of the evaluations recorded so far, if there are any. + + This is what an interrupted or crashed search leaves behind: the records + are already on disk, one per evaluation, and consolidating them costs + nothing next to the model runs they stand for. A failure to write the table + is reported and swallowed, since the exception that brought the search down + is the one the caller has to see. + """ + try: + records = _read_records(evaluations_dir) + if records: + _write_evaluations_csv(evaluations_csv, records) + logger.warning( + "The search ended early; the %d evaluation(s) it recorded are in %s.", + len(records), + evaluations_csv, + ) + except Exception as error: # pragma: no cover - the original failure must win + logger.warning("The evaluations of the interrupted search could not be written: %s", error) + + +def _compared_steps(*, first_step: int, last_step: int, spinup_steps: int) -> list[int]: + """Return the simulated time steps the efficiency is computed on, in order.""" + return [step for step in range(first_step, last_step + 1) if step > spinup_steps] + + +def _check_station_selection( + stations: tuple[str, ...] | None, + observed: Series, + sampled_ids: set[str] | None, +) -> None: + """Refuse a station selection that names none of the stations of the run. + + The selection is applied inside every evaluation, so a misspelt id would + make every one of them fail and the whole budget would be spent finding + out. The comparable stations are the ones the observed series and the run + have in common; under the ``zones`` aggregation the parent does not know + the ids of the run, so the selection is only checked against the + observations. + + :param stations: The ids the objective averages, or ``None``. + :type stations: tuple[str, ...] | None + + :param observed: The observed series, as the parent read it. + :type observed: rubem.calibration.objective.Series + + :param sampled_ids: The ids the run will sample, or ``None`` when they are + not known before the run. + :type sampled_ids: set[str] | None + + :raises CalibrationError: If none of the selected ids is comparable. + """ + if stations is None: + return + comparable = set(observed.stations) + if sampled_ids is not None: + comparable &= sampled_ids + unknown = [station for station in stations if station not in comparable] + if len(unknown) == len(stations): + raise CalibrationError( + f"The station(s) {', '.join(stations)} of the objective are none of the " + f"station(s) {', '.join(sorted(comparable)) or 'none'} the observed series and " + "the run have in common; every evaluation would fail." + ) + if unknown: + logger.warning( + "The objective names %d station(s) the comparison does not have (%s); it is " + "averaged over the remaining ones.", + len(unknown), + ", ".join(unknown), + ) + + +def _starting_point(space: DecisionSpace, parameters: dict[str, float]) -> np.ndarray: + """Return the ``x0`` of the search: the configuration, inside the bounds of the run. + + The starting point is the parameter set of the configuration being + calibrated. A bound narrowed for this run may exclude it, and SciPy refuses + a starting point outside the bounds, so every coordinate that falls outside + is moved onto the bound it crosses and the move is reported: the run then + starts from the admissible point closest to the configuration instead of + ending before its first evaluation. + + :param space: The decision vector of the run. + :type space: rubem.calibration.parameters.DecisionSpace + + :param parameters: The calibration parameters of the configuration. + :type parameters: dict[str, float] + + :return: The free values of the starting point, in the order of the space. + :rtype: numpy.ndarray + """ + x0 = space.from_parameters(parameters) + for index, (name, (minimum, maximum)) in enumerate( + zip(space.free_names, space.bounds, strict=True) + ): + value = float(x0[index]) + moved = min(max(value, minimum), maximum) + if moved != value: + x0[index] = moved + logger.warning( + "The configuration has %s = %g, outside the bound [%g, %g] of this " + "calibration; the search starts from %g instead.", + name, + value, + minimum, + maximum, + moved, + ) + if not space.is_admissible(x0): + logger.warning( + "The parameters of the configuration are not an admissible candidate of this " + "calibration; the search starts from the rest of the initial population." + ) + return x0 + + +def _report_decision_space(space: DecisionSpace, stations: tuple[str, ...] | None) -> None: + """Log what this run searches: the fixed parameters, the narrowed bounds and the stations.""" + if space.fixed: + logger.info( + "Fixed, not searched: %s.", + ", ".join(f"{name} = {value:g}" for name, value in space.fixed.items()), + ) + narrowed = [ + (name, bound) + for name, bound in zip(space.free_names, space.bounds, strict=True) + if bound != variable_range(name) + ] + if narrowed: + logger.info( + "Searched in a narrowed range: %s.", + ", ".join( + f"{name} in [{minimum:g}, {maximum:g}]" for name, (minimum, maximum) in narrowed + ), + ) + if stations is not None: + logger.info("The objective averages the station(s) %s.", ", ".join(stations)) + + +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 _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 _observed_window(observed: Series, compared_steps: list[int]) -> np.ndarray: + """Return the rows of the observed series that fall on the compared steps.""" + steps = np.asarray(observed.steps, dtype=np.int64) + return np.flatnonzero(np.isin(steps, np.asarray(compared_steps, dtype=np.int64))) + + +def _in_selection(station: str, stations: tuple[str, ...] | None) -> str: + """Return whether a station enters the objective, as the tables spell it.""" + return "true" if stations is None or station in stations else "false" + + +def _write_observed_csv( + path: Path, + observed: Series, + compared_steps: list[int], + stations: tuple[str, ...] | None, +) -> None: + """Write what the observations offer on the compared steps, one row per station. + + The row of a station carries how many of the compared steps it observes + (``pairs_in_window``), how many values on those steps the mask of + :func:`rubem.calibration.objective.valid_mask` rejected as gaps + (``dropped``, the negative values among them) and the summary statistics of + the values that were kept. It is written before the search: a station whose + gauge is nearly empty over the simulated window explains a mean efficiency + afterwards, and it is cheaper to read it first. + """ + window = _observed_window(observed, compared_steps) + with path.open("w", encoding="utf-8", newline="") as table: + writer = csv.writer(table, delimiter=STATION_SEPARATOR) + writer.writerow(OBSERVED_COLUMNS) + for station in sorted(observed.stations): + values = np.asarray(observed.stations[station], dtype=np.float64)[window] + mask = valid_mask(values) + kept = values[mask] + writer.writerow( + [ + station, + _in_selection(station, stations), + int(kept.size), + int(values.size - kept.size), + _cell(float(np.mean(kept)) if kept.size else None), + _cell(float(np.std(kept, ddof=1)) if kept.size > 1 else None), + _cell(float(np.min(kept)) if kept.size else None), + _cell(float(np.max(kept)) if kept.size else None), + ] + ) + + +def _report_observed_coverage( + observed: Series, + compared_steps: list[int], + stations: tuple[str, ...] | None, +) -> None: + """Log how much of the compared window the observations cover, per station.""" + window = _observed_window(observed, compared_steps) + covered = np.asarray(observed.steps, dtype=np.int64)[window] + if covered.size: + span = f"the step(s) {int(covered.min())} to {int(covered.max())}" + else: + span = "no step of it" + progress.info( + "The calibration compares %d time step(s); the observed series covers %d of them, %s.", + len(compared_steps), + int(covered.size), + span, + ) + for station in sorted(observed.stations): + values = np.asarray(observed.stations[station], dtype=np.float64)[window] + kept = int(np.count_nonzero(valid_mask(values))) + progress.info( + "Station %s%s: %d observed value(s) on the compared steps, %d dropped as gaps.", + station, + "" if stations is None or station in stations else " (not in the objective)", + kept, + int(values.size) - kept, + ) + + +def _write_stations_csv( + path: Path, + station_metrics: dict[str, Any], + stations: tuple[str, ...] | None, +) -> None: + """Write the goodness of fit of the best candidate, one row per station. + + The statistics are the ones its own evaluation recorded, not a + recomputation: the table and the record of that evaluation therefore always + agree. A station outside the selection is measured like any other and + marked as such, which is what makes the stations left out of the objective a + validation of the calibration. + """ + with path.open("w", encoding="utf-8", newline="") as table: + writer = csv.writer(table, delimiter=STATION_SEPARATOR) + writer.writerow(STATION_COLUMNS) + for station in sorted(station_metrics): + measured = station_metrics[station] or {} + writer.writerow( + [ + station, + _in_selection(station, stations), + int(measured.get("pairs") or 0), + # The remaining columns are the statistics, under the very + # names the record spells them with. + *(_cell(measured.get(name)) for name in STATION_COLUMNS[3:]), + ] + ) + + +def _write_best_series( + path: Path, + *, + document: dict, + base_dir: str | None, + parameters: dict[str, float], + variable: str, + observed: Series, + spinup_steps: int, + temp_dir: Path, + workers: int, + evaluations_csv: Path, +) -> None: + """Run the best candidate once more and write its series beside the observed one. + + The search keeps no output: every evaluation writes into a temporary + directory that is removed when it ends, which is what keeps a calibration of + thousands of runs from filling a disk. The candidate that won is therefore + run one last time, in a worker of its own, so that the series behind the + reported efficiency can be plotted against the observations. The steps and + the stations of the table are the ones the objective compared, from the same + alignment, and a value either series does not offer is left as an empty + cell. + + :raises CalibrationError: If the run of the best candidate failed. + """ + output_dir = Path(tempfile.mkdtemp(prefix="best-", dir=str(temp_dir))) + executor = _pool(workers) + try: + table = executor.submit( + simulate_best, document, base_dir, parameters, str(output_dir), variable + ).result() + simulated = read_series(table) + except Exception as error: + raise CalibrationError( + "The best candidate of the search could not be run again to keep its series: " + f"{type(error).__name__}: {error}. The calibration itself finished; its " + f"evaluations are in {evaluations_csv}." + ) from error + finally: + executor.shutdown(wait=True) + shutil.rmtree(output_dir, ignore_errors=True) + + alignment = align_series(simulated, observed, spinup_steps) + columns = [ + column + for station in alignment.stations + for column in (f"observed_{station}", f"simulated_{station}") + ] + with path.open("w", encoding="utf-8", newline="") as best: + writer = csv.writer(best, delimiter=STATION_SEPARATOR) + writer.writerow(["step", *columns]) + for row, step in enumerate(alignment.steps.tolist()): + cells: list[str] = [] + for station in alignment.stations: + cells.append( + _observation(observed.stations[station], alignment.observed_index[row]) + ) + cells.append( + _observation(simulated.stations[station], alignment.simulated_index[row]) + ) + writer.writerow([int(step), *cells]) + + +def _observation(values: np.ndarray, row: int) -> str: + """Return one value of a series, or an empty cell when the mask rejects it.""" + value = float(np.asarray(values, dtype=np.float64)[row]) + return _cell(value) if bool(valid_mask(np.asarray([value]))[0]) else "" + + +def _population_size(settings: CalibrationSettings, dimension: int = len(FREE_PARAMETERS)) -> 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 * dimension)`` members, with + ``dimension`` the number of free parameters, and ``sobol`` alone 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 * dimension) + 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 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) + 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", ""), + record.get("started_at", ""), + *(_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. + + A record without an objective is not a candidate for the best one: the + records are written by other processes, and the parent reads whatever is in + the directory rather than assuming every file there is complete. + """ + evaluated = [ + record + for record in records + if not record.get("error") and record.get("objective") is not None + ] + 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, + space: DecisionSpace, + workers: int, + temp_dir: Path, + population_size: int, + artifacts: dict[str, str], +) -> None: + """Write the summary of the calibration. + + The ``artifacts`` object names every file the run wrote next to this one, so + that the summary is the one file a reader has to open to find the rest. + """ + 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, space, workers, temp_dir), + "artifacts": artifacts, + } + path.write_text(json.dumps(document, indent=2) + "\n", encoding="utf-8") + + +def _settings_document( + settings: CalibrationSettings, space: DecisionSpace, 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, the temporary directory and the bounds 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: the bounds are + the effective ones, one entry per searched parameter, the ranges of the + application settings narrowed by whatever the caller asked for. + """ + 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, + "bounds": { + name: [minimum, maximum] + for name, (minimum, maximum) in zip(space.free_names, space.bounds, strict=True) + }, + "fixed": dict(space.fixed), + "stations": list(settings.stations) if settings.stations is not None else None, + } + + +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_runner.py b/tests/unit/calibration/test_runner.py new file mode 100644 index 00000000..856e65a5 --- /dev/null +++ b/tests/unit/calibration/test_runner.py @@ -0,0 +1,1233 @@ +"""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 concurrent.futures.process import BrokenProcessPool +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.objective import Series +from rubem.calibration.parameters import ( + CALIBRATION_PARAMETERS, + FREE_PARAMETERS, + bounds, + parameters_to_vector, +) +from rubem.calibration.runner import ( + EVALUATION_COLUMNS, + OBSERVED_COLUMNS, + STATION_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 ( + Aggregation, + 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)) + + +def read_station_table(path): + """Read one of the ``;``-separated per-station tables of a run directory.""" + with path.open(encoding="utf-8", newline="") as table: + return list(csv.DictReader(table, delimiter=";")) + + +@pytest.fixture +def progress_log(caplog): + """Capture the ``rubem.progress`` lines whatever the logging configuration is. + + The command line detaches that logger from the root one so that its lines + are printed verbatim, and a test that ran the command line before this one + may have left it detached; attaching the capture handler to the logger + itself is what makes the capture independent of that. + """ + captured = logging.getLogger("rubem.progress") + handlers, level, propagate = captured.handlers[:], captured.level, captured.propagate + captured.handlers = [caplog.handler] + captured.setLevel(logging.INFO) + captured.propagate = False + caplog.set_level(logging.INFO) + try: + yield caplog + finally: + captured.handlers[:] = handlers + captured.setLevel(level) + captured.propagate = propagate + + +@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_every_artifact(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.stations_csv.is_file() + assert result.best_series.is_file() + assert (result.run_dir / "observed.csv").is_file() + assert result.best_series.name == "best_arn.csv" + assert result.message + assert result.generations == 1 + + @pytest.mark.unit + def test_the_summary_names_every_artifact_it_wrote(self, calibrated): + _, result = calibrated + artifacts = json.loads(result.result_json.read_text(encoding="utf-8"))["artifacts"] + + assert Path(artifacts["evaluations_csv"]) == result.evaluations_csv + assert Path(artifacts["observed_csv"]) == result.run_dir / "observed.csv" + assert Path(artifacts["stations_csv"]) == result.stations_csv + assert Path(artifacts["best_series"]) == result.best_series + assert Path(artifacts["calibrated_config"]) == result.calibrated_config + + @pytest.mark.unit + def test_the_observations_are_summarized_before_the_search(self, calibrated): + _, result = calibrated + rows = read_station_table(result.run_dir / "observed.csv") + + assert list(rows[0]) == list(OBSERVED_COLUMNS) + assert [row["station"] for row in rows] == ["1", "2"] + for row in rows: + # The observed file is the table the configuration itself wrote: + # every compared step is observed and nothing is a gap. + assert row["in_selection"] == "true" + assert row["pairs_in_window"] == str(TIMESTEPS) + assert row["dropped"] == "0" + assert float(row["min"]) <= float(row["mean"]) <= float(row["max"]) + assert float(row["std"]) >= 0.0 + + @pytest.mark.unit + def test_the_stations_of_the_best_candidate_are_a_table_of_their_own(self, calibrated): + _, result = calibrated + rows = read_station_table(result.stations_csv) + + assert list(rows[0]) == list(STATION_COLUMNS) + assert [row["station"] for row in rows] == ["1", "2"] + for row in rows: + # The optimum reproduces the observations exactly. + assert row["in_selection"] == "true" + assert int(row["pairs"]) == TIMESTEPS + assert float(row["nse"]) == pytest.approx(1.0, abs=1e-12) + assert float(row["rmse"]) == pytest.approx(0.0, abs=1e-9) + assert float(row["r"]) == pytest.approx(1.0, abs=1e-9) + assert float(row["mean_observed"]) == pytest.approx(float(row["mean_simulated"])) + assert float(row["std_observed"]) == pytest.approx(float(row["std_simulated"])) + + @pytest.mark.unit + def test_the_series_of_the_best_candidate_is_kept_beside_the_observed_one(self, calibrated): + _, result = calibrated + rows = read_station_table(result.best_series) + + assert list(rows[0]) == [ + "step", + "observed_1", + "simulated_1", + "observed_2", + "simulated_2", + ] + assert [int(row["step"]) for row in rows] == list(range(1, TIMESTEPS + 1)) + for row in rows: + for station in ("1", "2"): + # The best candidate is the configuration that wrote the + # observations, so the two columns of a station are one series. + assert float(row[f"simulated_{station}"]) == pytest.approx( + float(row[f"observed_{station}"]), abs=1e-9 + ) + + @pytest.mark.unit + def test_the_table_of_the_best_candidate_is_read_back_by_the_reader(self, calibrated): + from rubem.calibration.objective import read_series + + _, result = calibrated + table = read_series(result.best_series) + + # It is written in the layout the readers accept, header included. + assert table.steps.tolist() == list(range(1, TIMESTEPS + 1)) + assert sorted(table.stations) == ["observed_1", "observed_2", "simulated_1", "simulated_2"] + + @pytest.mark.unit + def test_every_row_of_the_table_carries_the_moment_it_started(self, calibrated): + import datetime + + _, result = calibrated + rows = read_evaluations(result.evaluations_csv) + + assert "started_at" in EVALUATION_COLUMNS + moments = [datetime.datetime.fromisoformat(row["started_at"]) for row in rows] + assert all(moment.tzinfo is not None for moment in moments) + # The rows are ordered by the candidate they evaluated; the moments are + # what puts them back in the order the evaluations were made in. + assert sorted(moments) != [] and len(moments) == len(rows) + + +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, progress_log, capfd + ): + evaluations = tmp_path / "evaluations" + evaluations.mkdir() + (evaluations / "one.json").write_text("{}", encoding="utf8") + callback = _Progress(evaluations) + + stop = callback(scipy_optimize.OptimizeResult(fun=12.5, nit=3)) + + assert stop is False, "the callback never halts the search" + assert callback.generations == 3 + assert "Generation 3" in progress_log.text + assert "12.5" in progress_log.text + assert "1 evaluation" in progress_log.text + # A record without an objective is not the best one; a callback that + # read it as a candidate would end the search it is only watching. + assert "best NSE n/a" in progress_log.text + assert capfd.readouterr().out == "", "the library logs, it does not print" + + @pytest.mark.unit + def test_the_generation_line_carries_the_best_efficiency_of_the_records( + self, tmp_path, progress_log + ): + evaluations = tmp_path / "evaluations" + evaluations.mkdir() + for name, objective_value, nse in (("a", 4.0, 0.25), ("b", 1.0, 0.75)): + (evaluations / f"{name}.json").write_text( + json.dumps({"id": name, "objective": objective_value, "nse": nse, "error": None}), + encoding="utf8", + ) + + _Progress(evaluations)(scipy_optimize.OptimizeResult(fun=1.0, nit=2)) + + # SciPy reports the objective; the efficiency behind it is read from + # the records, which is where the workers write it. + assert "Generation 2" in progress_log.text + assert "best NSE 0.750000" in progress_log.text + assert "2 evaluation(s) recorded" in progress_log.text + + +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 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 + 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 + ): + # 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="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" + 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" + + +@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_a_negative_observation_is_counted_as_a_gap_before_the_search( + self, dataset, broken_search, progress_log + ): + rows = dataset.observed.read_text(encoding="utf8").splitlines() + first = rows[1].split(";") + # The gauge of station 1 has no reading at the first step, written as + # the negative marker an observation table marks a gap with. + first[1] = "-1.5" + gapped = dataset.observed.with_name("gapped.csv") + gapped.write_text("\n".join([rows[0], ";".join(first), *rows[2:]]) + "\n", encoding="utf8") + + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate(observed=gapped) + + summarized = { + row["station"]: row for row in read_station_table(dataset.run_dir / "observed.csv") + } + # The negative value never reaches the statistics, and both the table + # and the line the run prints count it as a gap of that station alone. + assert summarized["1"]["dropped"] == "1" + assert summarized["1"]["pairs_in_window"] == str(TIMESTEPS - 1) + assert float(summarized["1"]["min"]) > 0.0 + assert summarized["2"]["dropped"] == "0" + assert summarized["2"]["pairs_in_window"] == str(TIMESTEPS) + assert ( + f"Station 1: {TIMESTEPS - 1} observed value(s) on the compared steps, " + "1 dropped as gaps." in progress_log.text + ) + + @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): + 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 + def test_a_smaller_decision_space_needs_fewer_members(self): + # One fixed parameter takes a dimension out of the search, and the + # number of members follows it. + assert _population_size(CalibrationSettings(popsize=15), 7) == 128 + assert _population_size(CalibrationSettings(popsize=15, init="latinhypercube"), 7) == 105 + assert _population_size(CalibrationSettings(popsize=1, init="random"), 1) == 5 + + @pytest.mark.unit + @pytest.mark.parametrize("popsize", [1, 5, 15, 16]) + @pytest.mark.parametrize("init", ["sobol", "halton", "latinhypercube", "random"]) + @pytest.mark.parametrize("dimension", [7, 8]) + def test_the_announced_population_is_the_one_scipy_builds(self, popsize, init, dimension): + # 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: + # only ``sobol`` rounds the population up to a power of two, and the + # number of free parameters is a factor of it. + optimum = scipy_optimize.differential_evolution( + lambda vector: float(np.sum(vector)), + bounds()[:dimension], + init=init, + popsize=popsize, + maxiter=0, + polish=False, + rng=0, + ) + + assert len(optimum.population) == _population_size( + CalibrationSettings(popsize=popsize, init=init), dimension + ) + + @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" + + +@pytest.fixture(scope="class") +def calibrated_with_a_fixed_parameter(tmp_path_factory): + """One calibration in which ``x`` is pinned and only seven parameters are searched.""" + data = Dataset(tmp_path_factory.mktemp("calibration-fixed")) + free = np.delete(INIT, FREE_PARAMETERS.index("x"), axis=1) + return data, data.calibrate(fixed={"x": 0.0}, init=free, stations=("1", "2")) + + +class TestFixedParameters: + @pytest.mark.unit + def test_the_fixed_parameter_is_not_searched_and_is_in_every_candidate( + self, calibrated_with_a_fixed_parameter + ): + _, result = calibrated_with_a_fixed_parameter + rows = read_evaluations(result.evaluations_csv) + + assert set(result.best_parameters) == set(CALIBRATION_PARAMETERS) + assert result.best_parameters["x"] == 0.0 + assert rows, "the search ran" + assert [row["error"] for row in rows if row["error"]] == [] + # Every candidate the search proposed carries the pinned value, and the + # column is still in the table: the parameter left the decision vector, + # not the model. + assert {row["x"] for row in rows} == {"0.0"} + + @pytest.mark.unit + def test_the_summary_records_the_fixed_value_the_bounds_and_the_stations( + self, calibrated_with_a_fixed_parameter + ): + _, result = calibrated_with_a_fixed_parameter + settings = json.loads(result.result_json.read_text(encoding="utf-8"))["settings"] + + assert settings["fixed"] == {"x": 0.0} + assert settings["stations"] == ["1", "2"] + assert list(settings["bounds"]) == [name for name in FREE_PARAMETERS if name != "x"] + assert settings["bounds"]["rcd"] == [1.0, 10.0] + + @pytest.mark.unit + def test_the_calibrated_configuration_carries_the_fixed_value( + self, calibrated_with_a_fixed_parameter + ): + _, result = calibrated_with_a_fixed_parameter + configuration = ModelConfiguration(result.calibrated_config, validate_input=False) + + assert configuration.calibration_parameters.model_dump()["x"] == 0.0 + + @pytest.mark.unit + def test_a_fixed_weight_replaces_the_constraint_with_a_narrowed_bound( + self, dataset, broken_search + ): + free = np.delete(INIT, FREE_PARAMETERS.index("w_1"), axis=1) + + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate(fixed={"w_1": 0.7}, init=free) + + searched = [name for name in FREE_PARAMETERS if name != "w_1"] + # One free weight is a bound, not a constraint: the argument is left out + # instead of being handed over as ``None``. + assert "constraints" not in broken_search + assert len(broken_search["bounds"]) == len(searched) == 7 + assert broken_search["bounds"][searched.index("w_2")] == (0.0, pytest.approx(0.3)) + assert len(broken_search["x0"]) == 7 + + @pytest.mark.unit + def test_a_decision_space_the_settings_cannot_describe_is_refused(self, dataset, no_pool): + with pytest.raises(ValueError, match="cannot be fixed"): + dataset.calibrate(fixed={"x": 5.0}) + + assert not dataset.run_dir.exists() + assert no_pool == [] + + @pytest.mark.unit + def test_an_initial_population_of_the_wrong_width_is_refused(self, dataset, no_pool): + # ``INIT`` has one column per parameter of the whole search, and the + # fixed one has left the decision vector. + with pytest.raises(ValueError, match="one column per free parameter") as failure: + dataset.calibrate(fixed={"x": 0.0}) + + assert "(5, 8)" in str(failure.value) + assert "7 free parameter(s)" in str(failure.value) + assert not dataset.run_dir.exists() + assert no_pool == [] + + @pytest.mark.unit + def test_an_initial_population_that_is_not_a_table_of_members_is_refused( + self, dataset, no_pool + ): + with pytest.raises(ValueError, match="two-dimensional array"): + dataset.calibrate(init=INIT[0]) + + assert not dataset.run_dir.exists() + assert no_pool == [] + + +class TestNarrowedBounds: + @pytest.mark.unit + def test_the_search_is_given_the_narrowed_bounds(self, dataset, broken_search): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate(bounds={"rcd": (2.0, 5.0)}) + + expected = bounds() + expected[FREE_PARAMETERS.index("rcd")] = (2.0, 5.0) + assert broken_search["bounds"] == expected + + @pytest.mark.unit + def test_a_starting_point_outside_a_narrowed_bound_is_moved_onto_it( + self, dataset, broken_search, caplog + ): + with caplog.at_level(logging.WARNING, logger="rubem.calibration.runner"): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate(bounds={"alpha": (5.0, 8.0)}) + + # The configuration has alpha = 4.5, which this run does not search; + # SciPy refuses a starting point outside the bounds, so it is moved onto + # the bound and the move is reported instead of ending the run. + assert broken_search["x0"][FREE_PARAMETERS.index("alpha")] == 5.0 + assert "alpha" in caplog.text + assert "4.5" in caplog.text + + @pytest.mark.unit + def test_a_bound_that_does_not_narrow_the_settings_range_is_refused(self, dataset, no_pool): + with pytest.raises(ValueError, match="not a narrower range"): + dataset.calibrate(bounds={"rcd": (0.5, 12.0)}) + + assert not dataset.run_dir.exists() + assert no_pool == [] + + +class TestStationSelection: + @pytest.mark.unit + @pytest.mark.parametrize( + "stations", + [(), ("1", ""), ("1", 2), "12"], + ids=["empty", "a blank id", "not a string", "a bare string"], + ) + def test_a_selection_that_names_no_station_is_refused(self, dataset, no_pool, stations): + with pytest.raises(CalibrationError): + dataset.calibrate(stations=stations) + + assert not dataset.run_dir.exists() + assert no_pool == [] + + @pytest.mark.unit + def test_the_selected_stations_are_logged_before_the_search( + self, dataset, broken_search, caplog + ): + with caplog.at_level(logging.INFO, logger="rubem.calibration.runner"): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate(stations=("1",)) + + assert "the station(s) 1" in caplog.text + + +def _record_then_raise(error): + """A search that records one evaluation and then fails, without a worker. + + The records are what an interrupted calibration has to leave behind, and + the failure has to happen where the real one does, inside the search. The + evaluations directory is read out of the objective the search was handed, + which is where the parent put it. + """ + + def stub(function, search_bounds, **kwargs): + evaluations_dir = Path(function.keywords["context"].evaluations_dir) + record = { + "id": "deadbeef", + "pid": 4321, + "started_at": "2000-01-01T00:00:00+00:00", + "parameters": dict.fromkeys(CALIBRATION_PARAMETERS, 0.5), + "nse": 0.25, + "station_nse": {"1": 0.25}, + "station_metrics": {}, + "objective": 5625000.0, + "elapsed_seconds": 1.5, + "error": None, + } + (evaluations_dir / "deadbeef.json").write_text(json.dumps(record), encoding="utf8") + raise error + + return stub + + +class TestInterruptedSearch: + @pytest.mark.unit + def test_an_interrupted_search_still_writes_the_evaluations_it_made(self, dataset, monkeypatch): + monkeypatch.setattr( + scipy_optimize, "differential_evolution", _record_then_raise(KeyboardInterrupt()) + ) + + with pytest.raises(KeyboardInterrupt): + dataset.calibrate() + + # The evaluations are hours of model runs on a real basin; a + # calibration that has to be stopped is read from them like any other. + rows = read_evaluations(dataset.run_dir / "evaluations.csv") + assert [row["id"] for row in rows] == ["deadbeef"] + assert rows[0]["started_at"] == "2000-01-01T00:00:00+00:00" + assert rows[0]["nse"] == repr(0.25) + + @pytest.mark.unit + def test_a_crashed_search_still_writes_the_evaluations_it_made(self, dataset, monkeypatch): + monkeypatch.setattr( + scipy_optimize, "differential_evolution", _record_then_raise(RuntimeError("boom")) + ) + + with pytest.raises(RuntimeError, match="boom"): + dataset.calibrate() + + assert read_evaluations(dataset.run_dir / "evaluations.csv") + + @pytest.mark.unit + def test_a_search_that_recorded_nothing_leaves_no_table(self, dataset, broken_search): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate() + + assert not (dataset.run_dir / "evaluations.csv").exists() + + @pytest.mark.unit + def test_a_worker_killed_by_the_system_is_reported_with_the_advice(self, dataset, monkeypatch): + monkeypatch.setattr( + scipy_optimize, + "differential_evolution", + _record_then_raise(BrokenProcessPool("A process in the process pool was terminated")), + ) + + with pytest.raises(CalibrationError, match="--workers") as failure: + dataset.calibrate() + + message = str(failure.value) + assert "died before finishing" in message + assert "memory" in message + # The cause is kept, so the traceback still says what the pool reported. + assert isinstance(failure.value.__cause__, BrokenProcessPool) + assert read_evaluations(dataset.run_dir / "evaluations.csv") + + +class TestProgressReporting: + @pytest.mark.unit + def test_the_budget_the_coverage_the_generations_and_the_summary_are_reported( + self, dataset, progress_log + ): + result = dataset.calibrate() + + text = progress_log.text + # The budget, before anything is spent on it. + assert "8 free parameter(s)" in text + assert f"{len(INIT)} population member(s)" in text + assert "2 worker process(es)" in text + # The coverage of the compared window, per station. + assert f"compares {TIMESTEPS} time step(s)" in text + assert f"covers {TIMESTEPS} of them" in text + assert f"Station 1: {TIMESTEPS} observed value(s)" in text + assert "0 dropped as gaps" in text + # One line per generation, and the closing summary. + assert "Generation 1:" in text + assert "best NSE" in text + assert "Calibration finished after" in text + assert f"{result.evaluations} evaluation(s)" in text + + @pytest.mark.unit + def test_the_progress_is_logged_and_never_printed(self, dataset, progress_log, capfd): + dataset.calibrate() + + assert "Generation 1:" in progress_log.text + out, _ = capfd.readouterr() + assert out == "", "the library logs, the command line prints" + + +class TestStationsOfTheObjective: + @pytest.mark.unit + def test_a_station_left_out_of_the_objective_is_measured_and_marked(self, dataset): + result = dataset.calibrate(stations=("1",)) + + observed = { + row["station"]: row for row in read_station_table(result.run_dir / "observed.csv") + } + stations = {row["station"]: row for row in read_station_table(result.stations_csv)} + + assert observed["1"]["in_selection"] == "true" + assert observed["2"]["in_selection"] == "false" + # The station outside the objective is measured all the same: that is + # the validation half of a calibration/validation split. + assert stations["1"]["in_selection"] == "true" + assert stations["2"]["in_selection"] == "false" + assert float(stations["2"]["nse"]) == pytest.approx(1.0, abs=1e-12) + assert int(stations["2"]["pairs"]) == TIMESTEPS + assert json.loads(result.result_json.read_text(encoding="utf-8"))["settings"][ + "stations" + ] == ["1"] + + @pytest.mark.unit + def test_a_selection_of_stations_the_run_does_not_have_is_refused_before_the_search( + self, dataset, no_pool + ): + with pytest.raises(CalibrationError, match="none of the station") as failure: + dataset.calibrate(stations=("7", "8")) + + # A misspelt id would otherwise fail every evaluation of the run and + # spend the whole budget finding out. + assert "7, 8" in str(failure.value) + assert "1, 2" in str(failure.value) + assert not dataset.run_dir.exists() + assert no_pool == [] + + @pytest.mark.unit + def test_a_selection_that_names_one_unknown_station_warns_and_keeps_the_rest( + self, dataset, broken_search, caplog + ): + with caplog.at_level(logging.WARNING, logger="rubem.calibration.runner"): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate(stations=("1", "7")) + + assert "does not have (7)" in caplog.text diff --git a/tests/unit/calibration/test_worker.py b/tests/unit/calibration/test_worker.py new file mode 100644 index 00000000..35948e6c --- /dev/null +++ b/tests/unit/calibration/test_worker.py @@ -0,0 +1,383 @@ +"""One evaluation of a candidate: the model run, the record and the failure modes.""" + +import datetime +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, + simulate_best, +) +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, stations=None): + """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, + stations=stations, + ) + + 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()) == [] + + +class TestRecordedMetrics: + @pytest.mark.unit + def test_the_record_carries_the_moment_the_evaluation_started(self, dataset): + before = datetime.datetime.now(datetime.UTC) + + evaluate(dataset.vector, dataset.context()) + + (record,) = dataset.records() + started_at = datetime.datetime.fromisoformat(record["started_at"]) + assert started_at.tzinfo is not None, "the moment carries its time zone" + assert started_at.utcoffset() == datetime.timedelta(0) + # The wall-clock moment is what puts the records back in the order they + # were made in, so it has to be the moment of this evaluation. + assert before <= started_at <= datetime.datetime.now(datetime.UTC) + + @pytest.mark.unit + def test_the_record_carries_the_goodness_of_fit_of_every_station(self, dataset): + evaluate(dataset.vector, dataset.context()) + + (record,) = dataset.records() + metrics = record["station_metrics"] + + assert sorted(metrics) == ["1", "2"] + for measured in metrics.values(): + assert sorted(measured) == [ + "mean_observed", + "mean_simulated", + "nse", + "pairs", + "r", + "rmse", + "std_observed", + "std_simulated", + ] + # The observations are what this very configuration wrote, so the + # simulation reproduces them exactly. + assert measured["pairs"] == TIMESTEPS + assert measured["nse"] == pytest.approx(1.0, abs=1e-12) + assert measured["rmse"] == pytest.approx(0.0, abs=1e-9) + assert measured["mean_observed"] == pytest.approx(measured["mean_simulated"]) + + @pytest.mark.unit + def test_the_efficiencies_are_the_ones_of_the_statistics(self, dataset): + # ``station_nse`` stays for the readers that only want the efficiency; + # it is taken from the statistics, so the two can never disagree. + evaluate(dataset.vector, dataset.context()) + + (record,) = dataset.records() + + assert record["station_nse"] == { + station: measured["nse"] for station, measured in record["station_metrics"].items() + } + + @pytest.mark.unit + def test_a_rejected_candidate_still_carries_the_moment_and_no_statistics(self, dataset): + from rubem.calibration.parameters import FREE_PARAMETERS + + vector = dataset.vector.copy() + vector[FREE_PARAMETERS.index("w_1")] = 0.6 + vector[FREE_PARAMETERS.index("w_2")] = 0.6 + + evaluate(vector, dataset.context()) + + (record,) = dataset.records() + assert record["error"] == "inadmissible" + assert record["station_metrics"] == {} + assert datetime.datetime.fromisoformat(record["started_at"]).tzinfo is not None + + @pytest.mark.unit + def test_a_station_outside_the_objective_is_measured_all_the_same(self, dataset): + # The objective averages station 1 alone; station 2 is reported anyway, + # which is the validation half of a calibration/validation split. + value = evaluate(dataset.vector, dataset.context(stations=("1",))) + + assert value == pytest.approx(0.0, abs=1e-9) + (record,) = dataset.records() + assert sorted(record["station_metrics"]) == ["1", "2"] + assert record["nse"] == pytest.approx(1.0, abs=1e-12) + + @pytest.mark.unit + def test_a_selection_that_names_no_shared_station_fails_the_evaluation(self, dataset): + from rubem.calibration.objective import INADMISSIBLE_OBJECTIVE + + value = evaluate(dataset.vector, dataset.context(stations=("nowhere",))) + + assert value == INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert "none of the stations" in record["error"] + + +class TestSimulateBest: + @pytest.mark.unit + def test_it_runs_the_candidate_and_keeps_the_table_it_wrote(self, dataset, tmp_path): + from rubem.calibration.parameters import vector_to_parameters + + output_dir = tmp_path / "best" + output_dir.mkdir() + + table = simulate_best( + dataset.document, + dataset.configuration.base_dir, + vector_to_parameters(dataset.vector), + str(output_dir), + "arn", + ) + + # The directory is left in place, which is the whole point: the series + # of the chosen candidate has to be readable after the run. + assert Path(table).is_file() + assert Path(table).parent == output_dir + assert sorted(entry.name for entry in output_dir.iterdir()) == [ + "metadata.json", + "tss_arn.csv", + ] + series = read_series(table) + assert sorted(series.stations) == ["1", "2"] + assert series.steps.tolist() == list(range(1, TIMESTEPS + 1)) + + @pytest.mark.unit + def test_it_reproduces_the_series_of_the_configuration_it_is_given(self, dataset, tmp_path): + from rubem.calibration.parameters import vector_to_parameters + + output_dir = tmp_path / "best" + output_dir.mkdir() + + table = simulate_best( + dataset.document, + dataset.configuration.base_dir, + vector_to_parameters(dataset.vector), + str(output_dir), + "arn", + ) + + simulated = read_series(table) + observed = read_series(dataset.observed_path) + for station, values in observed.stations.items(): + assert simulated.stations[station].tolist() == pytest.approx(values.tolist()) + + @pytest.mark.unit + def test_a_variable_the_run_does_not_write_raises(self, dataset, tmp_path): + from rubem.calibration.parameters import vector_to_parameters + + output_dir = tmp_path / "best" + output_dir.mkdir() + + # Unlike an evaluation, this one has nothing to rank a failure behind. + with pytest.raises(Exception, match="nope"): + simulate_best( + dataset.document, + dataset.configuration.base_dir, + vector_to_parameters(dataset.vector), + str(output_dir), + "nope", + ) From 9750912a57b15d2836419cac5d34fb09209e41f8 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Fri, 18 Sep 2026 16:12:41 -0300 Subject: [PATCH 5/9] Add the rubem calibrate command The command takes the configuration, the observed series and the run directory, and exposes the search as options: the variable, the spin-up window, the generations, the population, the seed, the workers, the temporary directory, per-run bounds (--bound NAME=MIN:MAX), fixed parameters (--fix NAME=VALUE), the stations of the objective, the initialization, the strategy and the polishing. It prints the progress lines and a closing summary listing every artifact, and exits without a traceback on configuration and calibration errors. --- rubem/cli.py | 274 ++++++++++++- tests/unit/test_cli_calibrate.py | 677 +++++++++++++++++++++++++++++++ 2 files changed, 948 insertions(+), 3 deletions(-) create mode 100644 tests/unit/test_cli_calibrate.py diff --git a/rubem/cli.py b/rubem/cli.py index fb9eb5fe..eb402753 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``. @@ -10,6 +11,7 @@ import logging import logging.config import logging.handlers +import math import sys import warnings from collections.abc import Sequence @@ -21,7 +23,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 @@ -70,6 +72,122 @@ class SchemaFormat(StrEnum): legacy = "legacy" +class InitMethod(StrEnum): + """Initial populations the differential evolution can be started from.""" + + sobol = "sobol" + latinhypercube = "latinhypercube" + halton = "halton" + random = "random" + + +def _parse_bounds(values: list[str] | None) -> dict[str, tuple[float, float]]: + """Return the ``--bound NAME=MIN:MAX`` options as a mapping of ranges. + + Only the shape of the text is checked here; whether the range narrows the + range of the application settings is settled by the calibration itself, + which knows the parameters. + + :param values: The options as they were written, or ``None``. + :type values: list[str] | None + + :return: The ``(minimum, maximum)`` range of each named parameter. + :rtype: dict[str, tuple[float, float]] + + :raises typer.BadParameter: If an option is not ``NAME=MIN:MAX``, or if a + parameter is bounded twice. + """ + parsed: dict[str, tuple[float, float]] = {} + for entry in values or []: + name, separator, rest = entry.partition("=") + minimum, colon, maximum = rest.partition(":") + if not separator or not colon or not name.strip(): + raise typer.BadParameter( + f"'{entry}' is not a bound; write it as NAME=MIN:MAX, as in 'rcd=2:5'.", + param_hint="'--bound'", + ) + parsed[_unique(name.strip(), parsed, "--bound")] = ( + _as_number(minimum, entry, "--bound"), + _as_number(maximum, entry, "--bound"), + ) + return parsed + + +def _parse_fixed(values: list[str] | None) -> dict[str, float]: + """Return the ``--fix NAME=VALUE`` options as a mapping of pinned values. + + :param values: The options as they were written, or ``None``. + :type values: list[str] | None + + :return: The value each named parameter is pinned to. + :rtype: dict[str, float] + + :raises typer.BadParameter: If an option is not ``NAME=VALUE``, or if a + parameter is fixed twice. + """ + parsed: dict[str, float] = {} + for entry in values or []: + name, separator, value = entry.partition("=") + if not separator or not name.strip(): + raise typer.BadParameter( + f"'{entry}' is not a fixed parameter; write it as NAME=VALUE, as in 'x=0'.", + param_hint="'--fix'", + ) + parsed[_unique(name.strip(), parsed, "--fix")] = _as_number(value, entry, "--fix") + return parsed + + +def _parse_stations(value: str | None) -> tuple[str, ...]: + """Return the ``--stations ID[,ID...]`` option as the ids it names. + + :param value: The option as it was written, or ``None``. + :type value: str | None + + :return: The station ids, in the order they were given. + :rtype: tuple[str, ...] + + :raises typer.BadParameter: If the option names no station, or leaves one of + the comma-separated ids empty. + """ + if value is None: + return () + stations = tuple(station.strip() for station in value.split(",")) + if not all(stations): + raise typer.BadParameter( + f"'{value}' does not name the stations; write them as ID[,ID...], as in '1,2,3'.", + param_hint="'--stations'", + ) + return stations + + +def _as_number(text: str, entry: str, option: str) -> float: + """Return a number written on the command line. + + :raises typer.BadParameter: If the text is not a finite number. + """ + try: + value = float(text) + except ValueError as error: + raise typer.BadParameter( + f"'{text.strip()}' of '{entry}' is not a number.", param_hint=f"'{option}'" + ) from error + if not math.isfinite(value): + raise typer.BadParameter( + f"'{text.strip()}' of '{entry}' is not a finite number.", param_hint=f"'{option}'" + ) + return value + + +def _unique(name: str, parsed: dict, option: str) -> str: + """Return a parameter name that was not given twice on the command line. + + :raises typer.BadParameter: If the name is already in the mapping. + """ + if name in parsed: + raise typer.BadParameter(f"'{name}' is given more than once.", param_hint=f"'{option}'") + return name + + def _version_callback(value: bool) -> None: if value: print(f"RUBEM v{__release__}") @@ -178,6 +296,156 @@ 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, + bound: Annotated[ + list[str] | None, + typer.Option("--bound", help="Narrow the range of one parameter: NAME=MIN:MAX."), + ] = None, + fix: Annotated[ + list[str] | None, + typer.Option("--fix", help="Keep one parameter out of the search: NAME=VALUE."), + ] = None, + stations: Annotated[ + str | None, + typer.Option( + "--stations", + help="Station ids of the objective, comma-separated (default: every shared station).", + ), + ] = None, + init: Annotated[ + InitMethod, typer.Option("--init", help="Initial population of the search.") + ] = InitMethod.sobol, + strategy: Annotated[ + str, typer.Option("--strategy", help="Differential evolution strategy.") + ] = "best1exp", + polish: Annotated[ + bool, + typer.Option("--polish/--no-polish", help="Refine the best candidate with a local search."), + ] = False, +) -> None: + """Calibrate the model parameters against an observed series.""" + # The options that carry a small language of their own are read first, so + # that a typo costs nothing but the message naming the form expected. + bounds = _parse_bounds(bound) + fixed = _parse_fixed(fix) + selected_stations = _parse_stations(stations) + + 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, + "init": init.value, + "strategy": strategy, + "polish": polish, + } + if workers is not None: + settings_arguments["workers"] = workers + if temp_dir is not None: + settings_arguments["temp_dir"] = str(temp_dir) + if bounds: + settings_arguments["bounds"] = bounds + if fixed: + settings_arguments["fixed"] = fixed + if selected_stations: + settings_arguments["stations"] = selected_stations + + 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) + print(f"Observed summary: {result.run_dir / 'observed.csv'}", flush=True) + print(f"Stations table: {result.stations_csv}", flush=True) + print(f"Best candidate series: {result.best_series}", flush=True) + + logger.info("RUBEM successfully finished!") + + @config_app.command("schema") def config_schema( schema_format: Annotated[ @@ -258,7 +526,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..5a3f8c57 --- /dev/null +++ b/tests/unit/test_cli_calibrate.py @@ -0,0 +1,677 @@ +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", + best_series=Path(run_dir) / "best_arn.csv", + stations_csv=Path(run_dir) / "stations.csv", + ) + + 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 + assert settings.init == "sobol" + assert settings.strategy == "best1exp" + assert settings.polish is False + assert settings.bounds is None + assert settings.fixed is None + assert settings.stations 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), + f"Observed summary: {result.run_dir / 'observed.csv'}", + f"Stations table: {result.stations_csv}", + f"Best candidate series: {result.best_series}", + ): + 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", + "--bound", + "--fix", + "--stations", + "--init", + "--strategy", + "--polish", + "--no-polish", + ): + 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" + + +class TestCliCalibrateDecisionSpaceOptions: + @pytest.mark.unit + def test_the_bounds_the_fixed_parameters_and_the_stations_reach_the_settings( + self, config_path, observed_path, tmp_path, fake_calibration, restore_logging + ): + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + "--bound", + "rcd=2:5", + "--bound", + "alpha=1:10", + "--fix", + "x=0", + "--stations", + "1, 2 ,3", + "--init", + "latinhypercube", + "--strategy", + "rand1bin", + "--polish", + ] + ) + + settings = fake_calibration.calls[0][3] + assert settings.bounds == {"rcd": (2.0, 5.0), "alpha": (1.0, 10.0)} + assert settings.fixed == {"x": 0.0} + # The ids are read as strings and the spacing around the commas is not + # part of them. + assert settings.stations == ("1", "2", "3") + assert settings.init == "latinhypercube" + assert settings.strategy == "rand1bin" + assert settings.polish is True + + @pytest.mark.unit + def test_no_polish_is_the_default_and_can_be_written_out( + self, config_path, observed_path, tmp_path, fake_calibration, restore_logging + ): + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + "--no-polish", + ] + ) + + assert fake_calibration.calls[0][3].polish is False + + @pytest.mark.unit + @pytest.mark.parametrize( + ("arguments", "message"), + [ + (["--bound", "rcd=2"], "NAME=MIN:MAX"), + (["--bound", "rcd:2:5"], "NAME=MIN:MAX"), + (["--bound", "=2:5"], "NAME=MIN:MAX"), + (["--bound", "rcd=a:5"], "is not a number"), + (["--bound", "rcd=2:nan"], "is not a finite number"), + (["--bound", "rcd=2:5", "--bound", "rcd=3:4"], "more than once"), + (["--fix", "x"], "NAME=VALUE"), + (["--fix", "=0"], "NAME=VALUE"), + (["--fix", "x=zero"], "is not a number"), + (["--fix", "x=0", "--fix", "x=1"], "more than once"), + (["--stations", ",1"], "ID[,ID...]"), + (["--stations", ""], "ID[,ID...]"), + (["--init", "quasirandom"], "quasirandom"), + ], + ) + def test_a_malformed_option_is_rejected_with_the_form_it_expects( + self, + config_path, + observed_path, + tmp_path, + fake_calibration, + capsys, + restore_logging, + arguments, + message, + ): + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + *arguments, + ] + ) + + assert error.value.code == 2 + assert message in capsys.readouterr().err + assert not fake_calibration.calls + + +class TestCliCalibrateProgress: + @pytest.mark.unit + def test_the_progress_logger_is_routed_to_standard_output(self, capsys, restore_logging): + """The lines a calibration reports reach the terminal verbatim.""" + import logging + + from rubem.cli import setup_logging + + setup_logging() + logging.getLogger("rubem.progress").info("Generation 7: best objective 1.") + + captured = capsys.readouterr() + assert captured.out == "Generation 7: best objective 1.\n" + assert captured.err == "", "the progress is not a diagnostic" + + @pytest.mark.unit + @pytest.mark.slow + def test_a_real_calibration_prints_the_budget_and_one_line_per_generation( + self, tmp_path, capsys, restore_logging + ): + """The whole command, on the synthetic dataset, with a search of two generations. + + The observations are the series the configuration itself wrote, so the + search has an optimum to find and every evaluation succeeds. The + population is the smallest a Sobol' initialization builds for eight free + parameters, which keeps the run to a few dozen model runs. + """ + pytest.importorskip("scipy.optimize") + + from rubem.api import Model + + config = write_synthetic_dataset(str(tmp_path), timesteps=3) + config_file = tmp_path / "config.json" + config_file.write_text(json.dumps(config), encoding="utf8") + observed = tmp_path / "observed.csv" + observed.write_bytes(Model.from_config(config).run().time_series["arn"][0].read_bytes()) + run_dir = tmp_path / "calibration" + + main( + [ + "calibrate", + "-c", + str(config_file), + "--observed", + str(observed), + "-o", + str(run_dir), + "--maxiter", + "1", + "--popsize", + "1", + "--workers", + "2", + "--seed", + "1", + "--temp-dir", + str(tmp_path / "temp"), + ] + ) + + output = capsys.readouterr().out + assert "Calibrating 8 free parameter(s)" in output, output + assert "compares 3 time step(s)" in output, output + assert "Generation 1:" in output, output + assert "best NSE" in output, output + assert "Calibration finished after" in output, output + assert "Calibration finished successfully!" in output + assert (run_dir / "observed.csv").is_file() + assert (run_dir / "stations.csv").is_file() + assert (run_dir / "best_arn.csv").is_file() From c8db15394da76bcf947d975434e6c952aede4ed6 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Fri, 18 Sep 2026 16:12:41 -0300 Subject: [PATCH 6/9] Add the dataset marker and the calibration smoke test on the Ipojuca basin Tests marked dataset need the local copy of the published datasets, which RUBEM_DATASET_DIR names; when it is unset, or does not name a directory, they are skipped with the reason, so the suite stays self-contained elsewhere. The smoke test rebases the manifest of the Ipojuca basin onto the local copy, restricts the window to six monthly steps, generates the drainage network once and names it in the configuration (the accumulated runoff is routed on a network lddcreate does not reproduce between calls), takes the series of a base run as the observations, and checks that a one-generation differential evolution from five members finds the configuration again with an objective of exactly zero. --- pyproject.toml | 1 + tests/conftest.py | 60 +++++++ tests/dataset/__init__.py | 0 tests/dataset/test_calibrate_ipojuca.py | 200 ++++++++++++++++++++++++ 4 files changed, 261 insertions(+) create mode 100644 tests/dataset/__init__.py create mode 100644 tests/dataset/test_calibrate_ipojuca.py 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(): 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 f94913b09813148c8981ea34a53d42408e436eb8 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Fri, 18 Sep 2026 16:12:41 -0300 Subject: [PATCH 7/9] Document the calibration Add the calibration page: what the command fits, the objective with its masks and guards, the free parameters with their bounds, the per-run bounds and fixed parameters, the derived weight, the differential evolution settings and the evaluation budget, the observed series layouts and the fixed drainage network a calibration on the accumulated runoff needs, the artifacts, the process model, what a run prints and hints for large machines. Document the command in the user guide and name the Python entry point in the API page. --- doc/source/api.rst | 31 ++ doc/source/calibration.rst | 857 +++++++++++++++++++++++++++++++++++++ doc/source/index.rst | 1 + doc/source/userguide.rst | 94 +++- 4 files changed, 978 insertions(+), 5 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..ded9d882 --- /dev/null +++ b/doc/source/calibration.rst @@ -0,0 +1,857 @@ +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, what it reports while +it runs, 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 +at most 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, and ``--fix`` takes further parameters out +of the search for one run. + +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 **negative**, or when +it is at or above ``1e30``. + +Every quantity the two series carry is a flux or a storage and cannot be +negative, so a negative entry is not a measurement but a gap: the ``-9999`` the +model and the usual observation tables write, the ``-1`` some gauge records +use, or any other negative marker of their own. All of them are dropped alike, +and nothing has to be declared for it. Zero is a value and not a gap: a station +may well record no flow. + +The ``1e30`` 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. + +How many observed values the mask dropped is reported per station before the +search starts, on the terminal and in the ``dropped`` column of +:file:`observed.csv`, so that a gauge that is nearly empty over the simulated +window is seen before the model runs and not after them. + +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 **selected** stations that +do have one. Without ``--stations`` every shared station is selected. +``--stations 1,2,3`` names the ids the mean is taken over, and the stations it +leaves out are still sampled, still compared and still measured: they carry +``false`` in the ``in_selection`` column of :file:`observed.csv` and of +:file:`stations.csv`, and their goodness of fit is recorded by every +evaluation. That is the calibration and validation split — the search is +driven by some of the gauges, and the others say what the parameters it found +do at a gauge that had no say in them. + +A selection that names none of the comparable stations is refused before the +search starts, since every evaluation would fail on it; a selection that names +some unknown id keeps the ids that exist and reports the others once. + +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 +:file:`zones_mapping.csv`, so that check, and the check of the selected +stations, are made against the observed series alone. + +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:: + + 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 +`````````````````` + +``--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. ``--bound`` and ``--fix`` accept either spelling, as they accept + ``w1`` and ``w2`` for ``w_1`` and ``w_2``. + +The ninth parameter is derived rather than searched: + +.. math:: + :label: calibration-w3 + :nowrap: + + \[w_3 = 1 - \left(w_1 + w_2\right)\] + +so a search over both weights is 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. + +Narrowing a range for one run +````````````````````````````` + +``--bound NAME=MIN:MAX``, repeatable, searches one parameter in a range of the +caller's choosing instead of the range of the application settings: + +.. code-block:: console + + $ rubem calibrate -c project-config.json --observed observed.csv -o calibration \ + --bound alpha=1:10 --bound w_1=0.1:1 --bound rcd=1:5 + +A bound may only **narrow** the settings range, never widen it: its minimum +must lie at or above the minimum of the settings, its maximum at or below their +maximum, and the minimum below the maximum. Anything else is refused before the +search starts, because a candidate outside the settings range is a +configuration the loader would reject anyway. Bounding the derived ``w_3`` is +refused as well: bound ``w_1`` and ``w_2`` instead. So is bounding a parameter +that is fixed, which has no range left to search. + +Narrowing is what a known basin buys: the ranges of the settings are the ones +the model accepts at all, and a range that past work has already narrowed +concentrates the population where the answer is instead of spending members on +values that were never plausible. + +Keeping a parameter out of the search +````````````````````````````````````` + +``--fix NAME=VALUE``, repeatable, pins a parameter for the whole run: + +.. code-block:: console + + $ rubem calibrate -c project-config.json --observed observed.csv -o calibration --fix x=0 + +The coordinate **leaves the decision vector**: the search loses one dimension, +the optimizer never moves the parameter, and the population is no longer spent +exploring it. The pinned value is nonetheless one of the nine parameters every +candidate is run with, and it is in the best parameters, in the calibrated +configuration and under ``settings.fixed`` in :file:`result.json`. + +A fixed value must lie inside the range of the application settings, and at +least one parameter must be left free; both are checked before the search +starts. The fixed parameters, and any range a ``--bound`` narrowed, are logged +by the ``rubem.calibration`` logger when the run starts, so a logging +configuration that shows informational records replays exactly what was +searched. + +The two weights are the case where fixing changes more than the dimension: + +- **Both weights searched**: the pair is kept under the linear constraint + :math:`w_1 + w_2 \leq 1`. +- **One weight fixed**: there is no pair left to constrain, so the constraint + becomes a bound on the weight that is still free — its maximum drops to + :math:`1 - w_{fixed} - w_{3,min}`, which is what keeps the derived + :math:`w_3` inside its own range. If that leaves the free weight a single + value (``--fix w_1=1`` leaves :math:`w_2 = 0`), that weight is fixed at it + as well and leaves the decision vector, so the dimension and the budget + describe the search that runs; if it leaves no value at all, the run is + refused instead of searching an empty range. +- **Both weights fixed**: :math:`w_3` is a single number. It is derived once + and checked before the search; a pair that does not add up to a valid third + weight is refused. + +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. A coordinate of +that starting point which a narrowed bound excludes is moved onto the bound it +crosses, and the move is logged: the run starts from the admissible point +closest to the configuration rather than ending before its first evaluation. + +The search +---------- + +Settings +```````` + +The search is SciPy's ``differential_evolution``. Three of the settings that +decide how it moves are on the command line, because a run may have to +reproduce the settings of an earlier calibration; the budget, ``--maxiter`` and +``--popsize``, is `The evaluation budget`_ and ``--seed`` is below; the rest +describe the method and are fixed. + +``--strategy`` (default ``best1exp``) + How a trial vector is built. With the default, it comes 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. Every strategy name SciPy accepts is accepted here. + +``--init`` (default ``sobol``) + The initial population: ``sobol``, ``latinhypercube``, ``halton`` or + ``random``. The default is a Sobol' sequence over the bounds, which covers + the box far more evenly than a random sample, so the initial evaluations are + not wasted on clusters and gaps. It is also the one initialization that + resizes the population — see `The evaluation budget`_. + +``--polish`` / ``--no-polish`` (default ``--no-polish``) + Whether the best candidate is refined by a local search once the last + generation is over; see `Polishing`_. + +``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. + +``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`_. + +Polishing +````````` + +Polishing is **off by default**, and on request with ``--polish``. When it is +on, the best candidate of the last generation is handed to a local search +(L-BFGS-B, or trust-constr when both weights are searched and the linear +constraint is therefore present), which keeps the result only if it improves +the objective. + +The reason not to pay for it by default is the price of a gradient here. The +local search has no derivatives to work from, so it estimates them by finite +differences: one model run per free parameter, plus the point itself, for every +step it takes, and those runs are made one at a time, outside the batch a +generation is spread over. On the full search one such step is nine serial +simulations, where a generation is 128 parallel ones. What it estimates is a gradient of an +objective read back from a simulated table, which the model does not provide +reliably at that scale. The evaluations polishing makes are recorded like any +other, so they appear in :file:`evaluations.csv` and in ``nfev``. + +The evaluation budget +````````````````````` + +The number of population members is not ``--popsize``: it is + +.. math:: + :label: calibration-members + :nowrap: + + \[members = \max\left(5,\ popsize \times N\right)\] + +with :math:`N` the number of free parameters — eight, minus one for every +``--fix``. With ``--init sobol``, the default, that count is then rounded up to +the next power of two, because a Sobol' sequence is balanced only over a +power-of-two sample; the other three initializations use it as it is. 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 and no fixed parameter, ``--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. The same run with +``--init latinhypercube`` keeps its **120 members** and at most +:math:`120 \times 101 = 12\,120` 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. +Fixing parameters shrinks :math:`N` and with it the population: four fixed +parameters leave :math:`N = 4`, so ``--popsize 15`` gives +:math:`15 \times 4 = 60`, rounded up to **64 members**. + +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. Polishing works the other way and adds +evaluations after the last generation. + +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. That number, and not the +announced one, is what ran. + +Inputs +------ + +The observed series +``````````````````` + +``--observed`` takes one file, in either of two layouts, recognized from its +content: a first line carrying the ``;`` separator is read as the CSV table of +the model, anything else as a PCRaster time series. + +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 first cell of that header is the label of the step column: ``0`` in the +tables the model writes, or any non-numeric word in a table written by hand. A +first cell that is another number is a row of values, which means the file has +no header; it is refused, naming both layouts. + +The **PCRaster TSS layout** is its documented form: a title line such as +``timeseries scalar``, the number of columns, one line per column name +beginning with the time step column, and then one whitespace separated row per +time step: + +.. code-block:: console + + timeseries scalar + 4 + timestep + 1 + 2 + 3 + 1 12.4 3.1 0.8 + 2 15.9 4.0 1.2 + +The header is required. A time series file without it is refused rather than +read with columns numbered ``1``, ``2``, ...: the columns of a time series are +stations, nothing in a head-less file says which station a column belongs to, +and a numbering would quietly compare the observations of one gauge with the +simulation of another. Adding the header to a file that lacks one is a title, a +count and one line per column, and the refusal names the layout it expects. + +.. 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. + +Gaps are written as they come: ``-9999``, any other negative marker, an empty +cell or a word that is not a number. The mask of `Which values count`_ is what +decides, so nothing has to be converted before a file is handed over. + +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 +------------------ + +What the run directory, ``-o``/``--run-dir``, receives. The tables of the +stations use ``;`` as their separator, the one the model writes its own tables +with; :file:`evaluations.csv` is comma-separated. + +:file:`observed.csv` + Written **before** the search, one row per station of the observed series, + with the columns ``station``, ``in_selection``, ``pairs_in_window``, + ``dropped``, ``mean``, ``std``, ``min`` and ``max``. ``in_selection`` says + whether the station is one the objective averages over, which is every + station without ``--stations``; a station of the observed file the + configuration does not sample is listed here all the same, and never reaches + the objective. ``pairs_in_window`` is how many of the compared steps it + observes with a value the mask accepts, ``dropped`` how many values on those + steps the mask rejected as gaps, and the remaining columns summarize the + values that were kept. It is what a mean efficiency is + read against afterwards, and it costs nothing to look at first. + +:file:`evaluations.csv` + One row per evaluation, with the columns ``id``, ``pid``, ``started_at``, + ``alpha``, ``beta``, ``w_1``, ``w_2``, ``w_3``, ``rcd``, ``f``, + ``alpha_gw``, ``x``, ``nse``, ``objective``, ``elapsed_seconds`` and + ``error``. ``id`` identifies the evaluation, ``pid`` the process that made + it and ``started_at`` the wall-clock moment it began, in UTC; ``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 — and ``started_at`` is what puts them back in the order + they were made in, which is the order a convergence curve is drawn in. + +:file:`stations.csv` + The goodness of fit of the best candidate, one row per station the two + series share, with the columns ``station``, ``in_selection``, ``pairs``, + ``mean_observed``, ``std_observed``, ``mean_simulated``, ``std_simulated``, + ``r`` (Pearson correlation), ``rmse`` and ``nse``. The statistics are the + ones the evaluation of that candidate recorded, not a recomputation, so the + table and the record always agree. A station outside the selection is + measured like any other and marked as such: this is where the validation + gauges are read. + +:file:`best_.csv` + The series behind the reported efficiency: the column ``step``, then + ``observed_`` and ``simulated_`` for every shared station, one row + per compared step. The steps and the stations are the ones the objective + compared, and a value either series does not offer is left as an empty cell. + The search keeps no output of its own — every evaluation writes into a + temporary directory that is removed again — so the winning candidate is run + once more, after the search, to produce this table. It is the last thing the + run does before the summary, so a failure of that final run is reported with + :file:`stations.csv` and the calibrated configuration already written and + :file:`result.json` not written at all; the search itself has finished and + its table is complete. + +: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``, an ``artifacts`` object naming + every file above, and a ``settings`` object with ``variable``, + ``spinup_steps``, ``maxiter``, ``popsize``, ``seed``, ``workers``, + ``temp_dir``, ``init``, ``strategy``, ``mutation``, ``recombination``, + ``polish``, ``bounds``, ``fixed`` and ``stations``. The workers, the + temporary directory and the bounds are the values the run resolved and used + — the bounds one entry per searched parameter, the settings ranges already + narrowed — and 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_metrics``, the statistics of +:file:`stations.csv` for every station of that evaluation, and ``station_nse``, +the efficiencies alone — which is where to look when the mean hides one station +behaving differently from the others, for a candidate that is not the best one. + +.. 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. + +An interrupted run +`````````````````` + +A search that is interrupted, or that dies, after its first evaluation still +leaves the evaluations it made: the records are already on disk, and they are +consolidated into :file:`evaluations.csv` on the way out, with a message saying +where they are. Together with the :file:`observed.csv` written before the +search and the records themselves, that is the whole history of the run. What +needs the finished search is not written: :file:`result.json`, +:file:`stations.csv`, :file:`best_.csv` and the calibrated +configuration. Since a run directory is never reused, an interrupted run is +read from its table and a new directory is used for the next attempt. + +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. + +What a run prints +````````````````` + +The progress of a calibration goes to the ``rubem.progress`` logger, which the +command line prints on standard output, exactly as it does for a simulation; a +calibration started from Python says nothing until its host configures that +logger. Four kinds of line are written, in this order: + +.. code-block:: console + + The calibration compares 216 time step(s); the observed series covers 216 of them, the step(s) 13 to 228. + Station 1: 198 observed value(s) on the compared steps, 18 dropped as gaps. + Station 2 (not in the objective): 204 observed value(s) on the compared steps, 12 dropped as gaps. + 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). + Generation 1: best objective 1.60757e+06, best NSE 0.601240, 256 evaluation(s) recorded. + Generation 2: best objective 4.21033e+05, best NSE 0.794900, 384 evaluation(s) recorded. + Calibration finished after 2048 evaluation(s) in 16 generation(s): best objective 160757, best NSE 0.873210. + +The coverage lines say how much of the compared window the observations offer, +per station, and which stations the objective leaves out. The budget line is +the arithmetic of `The evaluation budget`_. A generation line reports the best +objective the optimizer holds and, beside it, the best efficiency among the +records the workers have written so far, with how many evaluations have been +recorded; a best objective of ``1e+30`` means no candidate has been evaluated +successfully yet. The closing line repeats the search as SciPy reports it. + +Everything else — the fixed parameters, the narrowed ranges, the stations of +the objective, the warnings about a starting point moved onto a bound, and the +summary of the failed evaluations — goes to the ``rubem.calibration`` logger. +The command line shows its warnings; its informational records need a logging +configuration that asks for them. + +When a worker is killed +``````````````````````` + +A worker that disappears before it finishes its evaluation ends the +calibration, and the most common reason is the system killing it for lack of +memory: every worker holds the rasters of one model run. The message says so, +advises a lower ``--workers`` and names the table of the evaluations made up to +that point, which is written on the way out. + +Definitions worth knowing +------------------------- + +Five things about the method that are easier to read once than to infer from a +table of results. + +The weights are derived exactly + ``w_3`` is computed as :math:`1 - (w_1 + w_2)`, a single subtraction, so the + three weights of a candidate add up to 1 exactly and no candidate is ever + run with weights that do not. There is no tolerance and no penalty term: a + pair of weights that would leave ``w_3`` outside its range is not scored at + all, it is rejected before the model runs. + +A time step is a month + The model computes a monthly water balance, so the steps the efficiency is + computed on are the months of ``SIM_TIME``, and a window of 228 steps is + nineteen years. The steps carry the numbers the model gives them, counted + from the alignment month of the raster series rather than from the start of + the simulated period, so the first compared step need not be ``1``. + ``--spinup-steps`` drops steps by that number, and the coverage line the run + prints says which steps are left. + +A ``point`` station should mark one cell + The value a time series carries for a station is the average over all the + cells that carry that id in the sample locations raster. With the ``point`` + aggregation a station is meant to be a gauge, so its id should mark exactly + one cell; an id painted over several cells is compared as the mean of them, + which is a different quantity from the discharge of a gauge. The + ``subcatchment`` and ``zones`` aggregations are the ones that average an + area on purpose. + +The reported best was really run + The best candidate is one the search evaluated: its row is in + :file:`evaluations.csv`, with the objective the run produced. Nothing is + extrapolated, interpolated or reconstructed from the population, and + :file:`stations.csv` repeats the statistics of that very evaluation. + +The efficiency of the winner, station by station + :file:`stations.csv` has it, one row per station, with ``in_selection`` + separating the gauges that drove the search from the ones that only watched. + The same numbers are in the JSON record of that evaluation, under + ``station_metrics`` and ``station_nse``, and :file:`best_.csv` + holds the two series the numbers come from, ready to plot. + +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, and a + job that asks for more workers than its memory allows is stopped by the + system killing one of them. + +``--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. + +Fix and narrow what is already known + ``--fix`` removes a dimension from the search and ``--bound`` shrinks one. + Both make a given budget go further on the parameters that are still open, + and both are recorded in :file:`result.json`, so a run calibrated with part + of its parameters pinned still says which ones and at what value. + +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`, :file:`observed.csv`, :file:`stations.csv`, + :file:`best_.csv` and :file:`result.json` are the record of how a + parameter set was obtained: the settings, the seed, the number of + evaluations, what the observations offered 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 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. + --bound Narrow the range of one parameter: + NAME=MIN:MAX. + --fix Keep one parameter out of the search: + NAME=VALUE. + --stations Station ids of the objective, comma- + separated (default: every shared station). + --init + Initial population of the search. [default: + sobol] + --strategy Differential evolution strategy. [default: + best1exp] + --polish / --no-polish Refine the best candidate with a local + search. [default: no-polish] + -h, --help Show this message and exit. + +``calibrate`` needs SciPy, which comes with ``pip install +"rubem[calibration]"``. It writes ``observed.csv`` (what the observations offer +at each station, before the search), ``evaluations.csv`` (one row per +evaluation), ``stations.csv`` (the goodness of fit of the best candidate at +each station), ``best_.csv`` (its series beside the observed one), +``result.json``, ``-calibrated.json`` and the ``evaluations/`` +directory of per-evaluation records into the ``--run-dir``, which must be empty +of earlier records. ``--bound`` and ``--fix`` are repeatable, one parameter +each; ``--stations`` names the gauges the objective averages, and the others +are still measured and reported. The method itself (the objective, the +parameter bounds, the derivation of ``w3``, how to budget the number of +evaluations, what the run prints and what it needs from the configuration) is +documented in :doc:`Calibration `. + +.. code-block:: console + + $ rubem calibrate -c project-config.json --observed observed.csv -o calibration \ + --spinup-steps 12 + Loading configuration and validating inputs... + Calibration started... + The calibration compares 216 time step(s); the observed series covers 216 of them, the step(s) 13 to 228. + Station 1: 198 observed value(s) on the compared steps, 18 dropped as gaps. + 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). + Generation 1: best objective 1.60757e+06, best NSE 0.601240, 256 evaluation(s) recorded. + Generation 16: best objective 160757, best NSE 0.873210, 2048 evaluation(s) recorded. + Calibration finished after 2048 evaluation(s) in 16 generation(s): best objective 160757, best NSE 0.873210. + 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 + Observed summary: /basins/ipojuca/calibration/observed.csv + Stations table: /basins/ipojuca/calibration/stations.csv + Best candidate series: /basins/ipojuca/calibration/best_arn.csv + Use ``-V`` or ``--version`` to get the version of the RUBEM. .. code-block:: console From e88ed5255b5de11cbc8919e921556ec04cedc542 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Fri, 18 Sep 2026 16:12:56 -0300 Subject: [PATCH 8/9] Record the calibration command in the changelog --- doc/source/changelog.rst | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index e39985a1..cebd8276 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -13,6 +13,24 @@ The format follows `Keep a Changelog `__. Added ````` +- Added ``rubem calibrate``: a differential evolution (SciPy) over the free + calibration parameters, with ``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 series given in the table + layout of the model's time series or as a PCRaster time series with its + header. Negative, ``-9999`` and PCRaster missing values are gaps; a spin-up + window, per-run bounds (``--bound``), fixed parameters (``--fix``), the + stations of the objective (``--stations``) and the initialization, strategy + and polishing of the search are options. Every candidate runs in a fresh + spawned worker through the Python API, and the run directory receives + ``observed.csv``, ``evaluations.csv``, ``stations.csv``, + ``best_.csv``, ``result.json`` and ``-calibrated.json``, + with progress lines on the terminal. SciPy comes with the optional extra + ``rubem[calibration]``. The calibration page of the documentation describes + the method, the evaluation budget, the fixed drainage network a calibration + on ``arn`` needs and the process model; a ``dataset`` pytest marker with the + ``RUBEM_DATASET_DIR`` convention carries a smoke test on the Ipojuca basin + (`#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 99023b4a48fdc4485a99c8e2c6ef06c9394486ef Mon Sep 17 00:00:00 2001 From: camila billerbeck <50674811+camilabillerbeck@users.noreply.github.com> Date: Tue, 22 Sep 2026 16:00:23 -0300 Subject: [PATCH 9/9] Add support for --allow-blocking-problems option in calibration - Introduced the --allow-blocking-problems flag in the CLI for the calibrate command to allow searches even with input validation errors. - Updated CalibrationSettings to include allow_blocking_problems parameter, and related documentation in calibration.rst, changelog.rst, and userguide.rst to reflect the new functionality. - Added unit tests to verify the behavior of the new option in both the calibration runner and CLI. --- doc/source/calibration.rst | 30 +++++++++++++++++++++ doc/source/changelog.rst | 6 ++++- doc/source/userguide.rst | 9 ++++++- rubem/calibration/runner.py | 20 +++++++++++++- rubem/cli.py | 11 ++++++++ tests/unit/calibration/test_runner.py | 38 +++++++++++++++++++++++++++ tests/unit/test_cli_calibrate.py | 20 ++++++++++++++ 7 files changed, 131 insertions(+), 3 deletions(-) diff --git a/doc/source/calibration.rst b/doc/source/calibration.rst index ded9d882..e97978d2 100644 --- a/doc/source/calibration.rst +++ b/doc/source/calibration.rst @@ -570,6 +570,36 @@ 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. +Inputs the validation rejects +````````````````````````````` + +The inputs are validated **once**, when the calibration loads the +configuration; the worker processes never revalidate them, since they do not +change during the search. A blocking problem therefore stops the calibration +before the first model run, which is the point: a search is hundreds or +thousands of runs, and an input the validation refuses would spend all of them +on a result that cannot be trusted. + +``--allow-blocking-problems`` searches anyway. The checks still run and every +problem is still reported — the non-blocking ones as warnings, the blocking +ones as errors — and the search then starts instead of stopping. It is the +counterpart of ``rubem run --allow-blocking-problems``, and there is no ``-s`` +on ``calibrate`` for it to conflict with. + +The value is written to :file:`result.json`, under ``settings``, because it +changes what the numbers beside it mean: the parameters were fitted on inputs +the validation rejected, and nothing else in the artifacts would say so. + +.. warning:: + + The option does not repair anything. A raster the validation blocks on still + reaches the model at every one of the evaluations, and the efficiency the + search maximizes is measured on what that produced. Whether the defect + reaches the compared stations is worth answering before the machine is + committed: a single ``rubem run --allow-blocking-problems`` writes the time + series of the stations, and values that are finite and plausible there say + the search has something real to fit. + Artifacts of a run ------------------ diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 5a425fcf..58d2a07a 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -40,7 +40,11 @@ Added not take it. ``ModelConfiguration(..., allow_blocking_problems=False)`` is the library counterpart and keeps the full list in ``ModelConfiguration.problems``; ``Model.from_file`` and ``Model.from_config`` of ``rubem.api`` take the same keyword and an isolated run - rebuilds the configuration with it + rebuilds the configuration with it. ``rubem calibrate --allow-blocking-problems``, and + ``CalibrationSettings(allow_blocking_problems=True)``, do the same for a calibration, + whose inputs are validated once before the search and never again by the workers; the + value is recorded under ``settings`` in :file:`result.json`, since the parameters it + reports were fitted on inputs the validation rejected (`#352 `__). - Added ``rubem.api``, the public Python surface: ``Model.from_file`` and ``Model.from_config`` load a configuration file or document, ``Model.run()`` diff --git a/doc/source/userguide.rst b/doc/source/userguide.rst index 0b0d7df2..37f2ebcf 100644 --- a/doc/source/userguide.rst +++ b/doc/source/userguide.rst @@ -1258,6 +1258,10 @@ Use ``-h`` or ``--help`` to get a brief description of each command and its argu best1exp] --polish / --no-polish Refine the best candidate with a local search. [default: no-polish] + --allow-blocking-problems Search even if the input validation finds + blocking problems; they are logged as + errors instead of stopping the + calibration. -h, --help Show this message and exit. ``calibrate`` needs SciPy, which comes with ``pip install @@ -1269,7 +1273,10 @@ each station), ``best_.csv`` (its series beside the observed one), directory of per-evaluation records into the ``--run-dir``, which must be empty of earlier records. ``--bound`` and ``--fix`` are repeatable, one parameter each; ``--stations`` names the gauges the objective averages, and the others -are still measured and reported. The method itself (the objective, the +are still measured and reported. ``--allow-blocking-problems`` is the +counterpart of the option of ``run``: the inputs are validated once, before the +search, and the option lets a configuration whose inputs carry blocking +problems be searched on anyway. The method itself (the objective, the parameter bounds, the derivation of ``w3``, how to budget the number of evaluations, what the run prints and what it needs from the configuration) is documented in :doc:`Calibration `. diff --git a/rubem/calibration/runner.py b/rubem/calibration/runner.py index 8ffc5501..e7610e5e 100644 --- a/rubem/calibration/runner.py +++ b/rubem/calibration/runner.py @@ -216,6 +216,15 @@ class CalibrationSettings: ``None``, the default, uses every station the observed series and the run have in common. :type stations: tuple[str, ...] | None + + :param allow_blocking_problems: Whether the calibration loads, and searches + on, a configuration whose inputs carry blocking problems. The checks + still run and every problem is still reported; the search then starts + instead of stopping. ``False``, the default, refuses such a + configuration, as a run of the model does. The value is recorded in + :file:`result.json`, since what the search fitted was measured on inputs + the validation rejected. + :type allow_blocking_problems: bool """ variable: str = "arn" @@ -233,6 +242,7 @@ class CalibrationSettings: bounds: dict[str, tuple[float, float]] | None = None fixed: dict[str, float] | None = None stations: tuple[str, ...] | None = None + allow_blocking_problems: bool = False @dataclass(frozen=True) @@ -439,7 +449,14 @@ def calibrate( ) config_file = as_path(config_path).absolute() - configuration = Model.from_file(config_file, validate_input=True).configuration + # The only validation of a whole calibration: the workers never revalidate, + # so this is where a blocking problem stops the search, and where + # ``allow_blocking_problems`` lets it start anyway. + configuration = Model.from_file( + config_file, + validate_input=True, + allow_blocking_problems=settings.allow_blocking_problems, + ).configuration file_v1 = _as_v1(configuration) _check_sample_locations(configuration, file_v1) @@ -1340,6 +1357,7 @@ def _settings_document( }, "fixed": dict(space.fixed), "stations": list(settings.stations) if settings.stations is not None else None, + "allow_blocking_problems": settings.allow_blocking_problems, } diff --git a/rubem/cli.py b/rubem/cli.py index d49716d7..cf19f25f 100644 --- a/rubem/cli.py +++ b/rubem/cli.py @@ -388,6 +388,16 @@ def calibrate( bool, typer.Option("--polish/--no-polish", help="Refine the best candidate with a local search."), ] = False, + allow_blocking_problems: Annotated[ + bool, + typer.Option( + "--allow-blocking-problems", + help=( + "Search even if the input validation finds blocking problems; " + "they are logged as errors instead of stopping the calibration." + ), + ), + ] = False, ) -> None: """Calibrate the model parameters against an observed series.""" # The options that carry a small language of their own are read first, so @@ -411,6 +421,7 @@ def calibrate( "init": init.value, "strategy": strategy, "polish": polish, + "allow_blocking_problems": allow_blocking_problems, } if workers is not None: settings_arguments["workers"] = workers diff --git a/tests/unit/calibration/test_runner.py b/tests/unit/calibration/test_runner.py index 856e65a5..b5f88f2a 100644 --- a/tests/unit/calibration/test_runner.py +++ b/tests/unit/calibration/test_runner.py @@ -213,6 +213,7 @@ def test_the_summary_carries_the_settings_and_the_population_size(self, calibrat assert summary["settings"]["polish"] is False assert summary["settings"]["temp_dir"] == str(data.temp_dir) assert summary["settings"]["init"] == INIT.tolist() + assert summary["settings"]["allow_blocking_problems"] is False @pytest.mark.unit def test_the_calibrated_configuration_loads_with_the_best_parameters(self, calibrated): @@ -776,6 +777,43 @@ def test_the_defaults_are_the_documented_ones(self): assert settings.recombination == 0.7 assert settings.polish is False assert settings.workers is None + assert settings.allow_blocking_problems is False + + @pytest.mark.unit + @pytest.mark.parametrize("allowed", [False, True]) + def test_allow_blocking_problems_reaches_the_configuration_loader( + self, dataset, monkeypatch, allowed + ): + """The single load of a calibration carries the option to ``Model.from_file``. + + The workers never revalidate, so this call is the only place the option + can have an effect, and the search itself is irrelevant to the question. + """ + seen = {} + original = Model.from_file + + def recording_from_file(path, **keywords): + seen.update(keywords) + return original(path, **keywords) + + def stop(*_args, **_keywords): + raise RuntimeError("the load is all this test exercises") + + monkeypatch.setattr(Model, "from_file", staticmethod(recording_from_file)) + # The search that follows the load is irrelevant here, so the next step + # of the calibration ends it. + monkeypatch.setattr(runner, "read_series", stop) + + with pytest.raises(RuntimeError, match="the load is all this test exercises"): + runner.calibrate( + dataset.config_file, + dataset.observed, + dataset.run_dir / "allowed", + CalibrationSettings(allow_blocking_problems=allowed), + ) + + assert seen["allow_blocking_problems"] is allowed + assert seen["validate_input"] is True @pytest.fixture(scope="class") diff --git a/tests/unit/test_cli_calibrate.py b/tests/unit/test_cli_calibrate.py index 5a3f8c57..984dcfd2 100644 --- a/tests/unit/test_cli_calibrate.py +++ b/tests/unit/test_cli_calibrate.py @@ -105,6 +105,7 @@ def test_the_paths_and_the_documented_defaults( assert settings.bounds is None assert settings.fixed is None assert settings.stations is None + assert settings.allow_blocking_problems is False @pytest.mark.unit def test_every_option_reaches_the_settings( @@ -556,6 +557,25 @@ def test_no_polish_is_the_default_and_can_be_written_out( assert fake_calibration.calls[0][3].polish is False + @pytest.mark.unit + def test_allow_blocking_problems_reaches_the_settings( + self, config_path, observed_path, tmp_path, fake_calibration, restore_logging + ): + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + "--allow-blocking-problems", + ] + ) + + assert fake_calibration.calls[0][3].allow_blocking_problems is True + @pytest.mark.unit @pytest.mark.parametrize( ("arguments", "message"),