From f5cc9e61fd3c7271d392aeaad4b396926fe4db34 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Thu, 17 Sep 2026 11:11:47 -0300 Subject: [PATCH 1/6] Add the decision space with per-run bounds and fixed parameters decision_space(fixed, bounds) builds the vector the search runs on: 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. The module-level helpers become thin wrappers over the default space. --- rubem/calibration/parameters.py | 491 ++++++++++++++++++---- tests/unit/calibration/test_parameters.py | 228 ++++++++++ 2 files changed, 637 insertions(+), 82 deletions(-) diff --git a/rubem/calibration/parameters.py b/rubem/calibration/parameters.py index de14d6ac..78c16efd 100644 --- a/rubem/calibration/parameters.py +++ b/rubem/calibration/parameters.py @@ -1,25 +1,42 @@ """The decision vector of the calibration and its translation to model parameters. -The differential evolution searches eight free parameters, ``alpha``, ``beta``, -``w_1``, ``w_2``, ``rcd``, ``f``, ``alpha_gw`` and ``x``; the slope factor -weight ``w_3`` is not searched but derived from the other two weights, since the -three weights of the weighted runoff coefficient must add up to 1. The search is -therefore eight-dimensional with the linear constraint ``w_1 + w_2 <= 1``, which -keeps the derived weight inside its own range. - -The bounds are the ones the application settings declare (``appsettings.json``, -section ``value_ranges.variables``), the same ranges +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. +configuration the model accepts. An override may only narrow them, never widen +them, for the same reason. -Nothing here imports SciPy at module level; :func:`weights_constraint`, the only -function that needs it, imports it when it is called. +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 typing import TYPE_CHECKING +from dataclasses import dataclass +from typing import TYPE_CHECKING, cast import numpy as np @@ -53,24 +70,387 @@ """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. + :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 [variable_range(name) for name in FREE_PARAMETERS] + 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 eight free values are read in the order of :data:`FREE_PARAMETERS` and - the slope factor weight is derived as ``w_3 = 1 - w_1 - w_2``, so that the - three weights always add up to 1. + 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 @@ -81,34 +461,14 @@ def vector_to_parameters(vector: Sequence[float] | np.ndarray) -> dict[str, floa :raises ValueError: If the vector does not have eight entries. """ - values = np.asarray(vector, dtype=np.float64).ravel() - if values.size != len(FREE_PARAMETERS): - raise ValueError( - f"The decision vector must have {len(FREE_PARAMETERS)} entries " - f"({', '.join(FREE_PARAMETERS)}), got {values.size}." - ) - parameters = {name: float(value) for name, value in zip(FREE_PARAMETERS, values, strict=True)} - # The two weights are added first and the sum is subtracted once. Taking - # them away one at a time, ``1 - w_1 - w_2``, rounds twice: for 20 of the - # 101 pairs of two-decimal weights that add up to 1, among them - # ``w_1 = 0.33, w_2 = 0.67``, it lands one ulp below zero, and - # :class:`rubem.configuration.calibration_parameters.CalibrationParameters` - # then refuses the run. With a single subtraction the weight is exactly - # zero at the boundary and never negative while ``w_1 + w_2 <= 1``. - parameters[DERIVED_PARAMETER] = 1.0 - (parameters["w_1"] + parameters["w_2"]) - return {name: parameters[name] for name in CALIBRATION_PARAMETERS} + 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. - This is how the configuration under calibration becomes the ``x0`` of the - differential evolution. Both spellings the configuration files use are - accepted, the canonical ``beta``, ``w_1``, ``w_2`` and the legacy ``b``, - ``w1``, ``w2``; the derived ``w_3`` is ignored, since it is not searched. - The straightforward way to obtain the mapping is to dump the loaded - :class:`rubem.configuration.calibration_parameters.CalibrationParameters` - of a configuration, which has already normalised the spellings. + 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] @@ -118,27 +478,14 @@ def parameters_to_vector(parameters: Mapping[str, float]) -> np.ndarray: :raises KeyError: If a free parameter is missing from the mapping. """ - normalised = {_ALIASES.get(name, name): value for name, value in parameters.items()} - missing = [name for name in FREE_PARAMETERS if name not in normalised] - if missing: - raise KeyError( - f"The calibration parameters are missing {', '.join(missing)}; " - f"the decision vector needs {', '.join(FREE_PARAMETERS)}." - ) - return np.asarray([float(normalised[name]) for name in FREE_PARAMETERS], dtype=np.float64) + 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. - A vector is admissible when every free value lies inside its bound, when - the two searched weights add up to at most 1 and when the derived ``w_3`` - lies inside its own bound. The derived weight is taken from - :func:`vector_to_parameters`, the very value the run would be given, so - that the guard can never disagree with the derivation and let through a - candidate - :class:`rubem.configuration.calibration_parameters.CalibrationParameters` - would then refuse. + 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 @@ -146,39 +493,19 @@ def is_admissible(vector: Sequence[float] | np.ndarray) -> bool: :return: ``True`` when the model can be run with this candidate. :rtype: bool """ - values = np.asarray(vector, dtype=np.float64).ravel() - if values.size != len(FREE_PARAMETERS): - return False - if not np.all(np.isfinite(values)): - return False - for value, (minimum, maximum) in zip(values, bounds(), strict=True): - if not minimum <= value <= maximum: - return False - - weights = float(values[FREE_PARAMETERS.index("w_1")] + values[FREE_PARAMETERS.index("w_2")]) - if weights > 1.0: - return False - derived = vector_to_parameters(values)[DERIVED_PARAMETER] - minimum, maximum = variable_range(DERIVED_PARAMETER) - return minimum <= derived <= maximum + return decision_space().is_admissible(vector) def weights_constraint() -> "LinearConstraint": """Return the linear constraint ``w_1 + w_2 <= 1`` on the decision vector. - The differential evolution receives it so that it never spends a model run - on a candidate whose slope factor weight would come out negative: SciPy - evaluates the constraint before the objective and skips the members that - violate it. + 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. """ - from scipy.optimize import LinearConstraint - - coefficients = np.zeros((1, len(FREE_PARAMETERS)), dtype=np.float64) - coefficients[0, FREE_PARAMETERS.index("w_1")] = 1.0 - coefficients[0, FREE_PARAMETERS.index("w_2")] = 1.0 - return LinearConstraint(coefficients, -np.inf, 1.0) + # 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 index 90d0db87..a2c209a7 100644 --- a/tests/unit/calibration/test_parameters.py +++ b/tests/unit/calibration/test_parameters.py @@ -8,6 +8,7 @@ DERIVED_PARAMETER, FREE_PARAMETERS, bounds, + decision_space, is_admissible, parameters_to_vector, vector_to_parameters, @@ -225,3 +226,230 @@ def test_the_constraint_evaluates_to_the_sum_of_the_weights(self): 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 0022bfea053d21557a5ae97958deacfc3467d048 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Thu, 17 Sep 2026 11:11:47 -0300 Subject: [PATCH 2/6] Treat negative observations as gaps, require headers and measure every station Every model output is a non-negative flux or storage, so a negative observed value is a gap. A ';' table without the "0;;..." header and a PCRaster time series without its header are refused naming the two accepted layouts, instead of numbering columns by position; a byte order mark is ignored. StationMetrics gives, per shared station, the valid pairs, the means and standard deviations of both series, the Pearson correlation, the RMSE and the NSE, and evaluate_series averages the NSE over a chosen selection of stations while measuring the others. --- rubem/calibration/objective.py | 369 +++++++++++++++++++---- tests/unit/calibration/test_objective.py | 323 +++++++++++++++++++- 2 files changed, 623 insertions(+), 69 deletions(-) diff --git a/rubem/calibration/objective.py b/rubem/calibration/objective.py index 7839851e..2272fb49 100644 --- a/rubem/calibration/objective.py +++ b/rubem/calibration/objective.py @@ -1,16 +1,17 @@ """The calibration objective: Nash-Sutcliffe efficiency over the sample stations. The module reads the observed and the simulated time series, masks the values -that stand for a missing observation, computes the Nash-Sutcliffe efficiency -(NSE) of each station and turns the mean efficiency into the objective the -optimizer minimizes. +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 and the objective are evaluated with -NumPy alone. +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 @@ -21,7 +22,12 @@ logger = logging.getLogger(__name__) MISSING_VALUES = (-9999.0,) -"""Values written by the model and by the usual observation tables for a gap.""" +"""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. @@ -49,9 +55,15 @@ def valid_mask(values: np.ndarray) -> np.ndarray: """Return the mask of the values that count as an observation. - A value counts when it is finite, is not one of :data:`MISSING_VALUES` and - is below ``1e30``, which covers the PCRaster missing value in any of its - round-trips. + 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 @@ -60,7 +72,7 @@ def valid_mask(values: np.ndarray) -> np.ndarray: :rtype: numpy.ndarray """ array = np.asarray(values, dtype=np.float64) - mask = np.isfinite(array) & (array < _MISSING_THRESHOLD) + mask = np.isfinite(array) & (array < _MISSING_THRESHOLD) & (array >= 0.0) for missing in MISSING_VALUES: mask &= array != missing return mask @@ -118,6 +130,121 @@ def nash_sutcliffe(simulated: np.ndarray, observed: np.ndarray) -> float | None: 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``. @@ -200,13 +327,36 @@ def _build(ids: list[str], steps: list[int], rows: list[list[str]], path: Path) ) +_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="utf8", newline="") as table: + 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) ] @@ -214,6 +364,7 @@ def _read_csv(path: Path) -> Series: 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]] = [] @@ -229,36 +380,61 @@ def _read_csv(path: Path) -> Series: 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. The model writes its own ``.tss`` files - without that header (see ``rubem.file._file_conversions.tss2csv``, which - reads every line of them as data); such a file is read as well, and its - stations are then numbered ``1``, ``2``, ... in column order, which is the - numbering of the sample locations of a usual dataset. + 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="utf8").splitlines() if line.strip()] + 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.") - ids: list[str] | None = None - data_start = 0 - if len(lines) > 1 and lines[1].strip().isdigit(): - columns = int(lines[1].strip()) - names = [line.strip() for line in lines[2 : 2 + columns]] - if len(names) == columns and columns >= 2: - ids = names[1:] - data_start = 2 + columns + 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 ids is None: - ids = [str(column) for column in range(1, len(cells))] if len(cells) != len(ids) + 1: raise ValueError( f"Row {number} of the time series file {path} has {len(cells)} column(s), " @@ -266,7 +442,7 @@ def _read_tss(path: Path) -> Series: ) steps.append(_as_step(cells[0], path, number)) rows.append(cells[1:]) - return _build(ids or [], steps, rows, path) + return _build(ids, steps, rows, path) def read_series(path: PathInput) -> Series: @@ -274,8 +450,14 @@ def read_series(path: PathInput) -> Series: The layout is recognised from the content: a first line carrying the ``;`` separator is the CSV table the model writes, anything else is read as a - PCRaster time series. Missing entries are kept as they are written; the - mask of :func:`valid_mask` is what decides whether a value counts. + 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] @@ -284,29 +466,49 @@ def read_series(path: PathInput) -> Series: :rtype: Series :raises FileNotFoundError: If the file does not exist. - :raises ValueError: If the file is empty, has no station column, repeats a - station id or has a row that does not match its header. + :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="utf8") + 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) -def evaluate_series( - simulated: Series, - observed: Series, - spinup_steps: int = 0, -) -> tuple[float, dict[str, float | None]]: - """Return the mean efficiency of a run and the efficiency of each station. +@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. The mean is - taken over the stations that have an efficiency, so that a station of zero - variance, reported as ``None``, does not drag the mean. + 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. @@ -320,11 +522,10 @@ def evaluate_series( :param spinup_steps: Number of initial steps to exclude, defaults to ``0``. :type spinup_steps: int, optional - :return: The mean efficiency and the efficiency of every shared station. - :rtype: tuple[float, dict[str, float | None]] + :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, share no time - step, or if no shared station has an efficiency. + :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: @@ -351,22 +552,76 @@ def evaluate_series( f"{observed_steps[observed_kept].tolist()}." ) - simulated_index = np.flatnonzero(simulated_kept)[simulated_rows] - observed_index = np.flatnonzero(observed_kept)[observed_rows] + 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 - station_nse: dict[str, float | None] = { - station: nash_sutcliffe( - simulated.stations[station][simulated_index], - observed.stations[station][observed_index], + :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 shared_ids + for station in alignment.stations } - values = [nse for nse in station_nse.values() if nse is not None] + 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 shared station " - f"({', '.join(shared_ids)}) has fewer than two valid pairs or observations " + "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)), station_nse + return float(np.mean(values)), metrics diff --git a/tests/unit/calibration/test_objective.py b/tests/unit/calibration/test_objective.py index bc76f666..92ebc5f2 100644 --- a/tests/unit/calibration/test_objective.py +++ b/tests/unit/calibration/test_objective.py @@ -1,5 +1,7 @@ """The calibration objective: masks, Nash-Sutcliffe efficiency, tables and alignment.""" +import dataclasses + import numpy as np import pytest @@ -8,10 +10,12 @@ MISSING_VALUES, PCRASTER_MISSING, Series, + StationMetrics, evaluate_series, nash_sutcliffe, objective, read_series, + station_metrics, valid_mask, ) @@ -20,7 +24,12 @@ FLOAT32_MISSING = float(np.float32(1e31)) # Every gap marker, used wherever a test drops a pair through one of them. -GAPS = (float("nan"), -9999.0, 1e31, FLOAT32_MISSING) +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): @@ -36,7 +45,19 @@ def series(steps, **stations): 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 + 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) @@ -251,17 +272,130 @@ def test_a_handwritten_pcraster_time_series(self, tmp_path): assert table.stations["22"].tolist() == [2.5, 4.5] @pytest.mark.unit - def test_a_pcraster_time_series_without_the_header(self, tmp_path): - # The form the model itself writes: every line is a data row and the - # stations are numbered in column order. + 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] - assert sorted(table.stations) == ["1", "2"] - assert table.stations["2"].tolist() == [2.5, 4.5] + + @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): @@ -306,7 +440,7 @@ def test_the_mean_over_the_stations(self): mean, per_station = evaluate_series(simulated, observed) - assert per_station == {"A": 1.0, "B": -3.0} + assert nse_of(per_station) == {"A": 1.0, "B": -3.0} assert mean == -1.0 @pytest.mark.unit @@ -318,7 +452,8 @@ def test_only_the_shared_steps_are_compared(self): mean, per_station = evaluate_series(simulated, observed) - assert per_station == {"A": 1.0} + 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 @@ -332,7 +467,7 @@ def test_the_spin_up_steps_are_excluded(self): mean, per_station = evaluate_series(simulated, observed, spinup_steps=2) - assert per_station == {"A": -3.0} + assert nse_of(per_station) == {"A": -3.0} assert mean == -3.0 @pytest.mark.unit @@ -342,7 +477,7 @@ def test_a_station_present_on_one_side_only_is_ignored(self): mean, per_station = evaluate_series(simulated, observed) - assert per_station == {"A": 1.0} + assert nse_of(per_station) == {"A": 1.0} assert mean == 1.0 @pytest.mark.unit @@ -354,7 +489,8 @@ def test_a_station_without_variance_is_reported_and_left_out_of_the_mean(self): mean, per_station = evaluate_series(simulated, observed) - assert per_station == {"A": 1.0, "B": None} + 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 @@ -394,3 +530,166 @@ def test_no_station_with_an_efficiency_is_an_error(self): 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 3239d65a0ebd450f513ce6897c1932dc8bf16627 Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Thu, 17 Sep 2026 11:11:47 -0300 Subject: [PATCH 3/6] Report the observations, the winner's series and the progress of a calibration The run directory receives observed.csv before the search (per station the selection, the pairs in the compared window, the dropped values and the summary statistics) and, after it, stations.csv with the winner's per-station metrics and best_.csv with its simulated series beside the observed one, from one more run of the best candidate in the pool. Every record carries started_at and station_metrics, the table the started_at column. The budget line, one line per generation with the best objective and NSE, and the final summary go through the rubem.progress logger the command line prints. An interrupted or failed search still consolidates evaluations.csv, and a worker killed by the system is reported with the advice to lower the number of workers. The decision space travels to the workers in the evaluation context. --- rubem/calibration/_worker.py | 153 +++++- rubem/calibration/runner.py | 714 +++++++++++++++++++++++--- tests/unit/calibration/test_runner.py | 532 ++++++++++++++++++- tests/unit/calibration/test_worker.py | 167 +++++- 4 files changed, 1449 insertions(+), 117 deletions(-) diff --git a/rubem/calibration/_worker.py b/rubem/calibration/_worker.py index edf32e89..b0d1e8b9 100644 --- a/rubem/calibration/_worker.py +++ b/rubem/calibration/_worker.py @@ -13,13 +13,17 @@ :meth:`rubem.api.Model.run_isolated`: the worker is already the isolation layer, and an isolated run inside it would nest one process pool in another. -Nothing here raises: a candidate outside the admissible region is rejected -without a run, and any failure of the run or of the evaluation is recorded and -turned into :data:`rubem.calibration.objective.INADMISSIBLE_OBJECTIVE`, so that -one failing candidate does not end the calibration. +: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 @@ -28,7 +32,8 @@ import time import uuid from collections.abc import Sequence -from dataclasses import dataclass +from dataclasses import dataclass, field +from datetime import UTC, datetime from pathlib import Path import numpy as np @@ -37,11 +42,12 @@ from .objective import ( INADMISSIBLE_OBJECTIVE, Series, + StationMetrics, evaluate_series, objective, read_series, ) -from .parameters import is_admissible, vector_to_parameters +from .parameters import DecisionSpace, decision_space logger = logging.getLogger(__name__) @@ -84,6 +90,16 @@ class EvaluationContext: 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 @@ -94,6 +110,8 @@ class EvaluationContext: 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: @@ -106,8 +124,9 @@ def evaluate(vector: Sequence[float] | np.ndarray, context: EvaluationContext) - objective. Every evaluation, successful or not, leaves one JSON record in the evaluations directory of the context. - :param vector: The candidate, the eight free values in the order of - :data:`rubem.calibration.parameters.FREE_PARAMETERS`. + :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. @@ -118,17 +137,21 @@ def evaluate(vector: Sequence[float] | np.ndarray, context: EvaluationContext) - candidate was rejected or its evaluation failed. :rtype: float """ + started_at = datetime.now(UTC).isoformat() started = time.perf_counter() try: - parameters = vector_to_parameters(vector) + parameters = context.space.to_parameters(vector) except ValueError as error: - _write_record(context, {}, None, {}, INADMISSIBLE_OBJECTIVE, 0.0, _describe(error)) + _write_record( + context, started_at, {}, None, {}, INADMISSIBLE_OBJECTIVE, 0.0, _describe(error) + ) return INADMISSIBLE_OBJECTIVE - if not is_admissible(vector): + if not context.space.is_admissible(vector): logger.debug("Rejecting an inadmissible candidate without running the model.") _write_record( context, + started_at, parameters, None, {}, @@ -140,13 +163,14 @@ def evaluate(vector: Sequence[float] | np.ndarray, context: EvaluationContext) - output_dir = tempfile.mkdtemp(prefix="evaluation-", dir=context.temp_dir) try: - nse, station_nse = _run_and_evaluate(context, parameters, output_dir) + 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, {}, @@ -159,16 +183,72 @@ def evaluate(vector: Sequence[float] | np.ndarray, context: EvaluationContext) - shutil.rmtree(output_dir, ignore_errors=True) value = objective(nse) - _write_record(context, parameters, nse, station_nse, value, time.perf_counter() - started, None) + _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, float | None]]: - """Run the candidate and return its mean efficiency and the efficiency per station.""" +) -> 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) @@ -176,7 +256,7 @@ def _run_and_evaluate( 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) + return evaluate_series(simulated, context.observed, context.spinup_steps, context.stations) def _time_series_table(result, variable: str) -> Path: @@ -198,18 +278,28 @@ 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 - temporary directory of the evaluation, every raster series is disabled (an - evaluation reads no raster, and writing them would dominate its cost) and - the time series are reduced to the calibrated variable, as CSV. The - aggregation the user configured is kept, since it decides which areas the - observed stations correspond to. + 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(context.document) + document = copy.deepcopy(base_document) document["model_calibration_parameters"] = { "alpha": parameters["alpha"], "b": parameters["beta"], @@ -231,7 +321,7 @@ def _derived_document( } output["time_series_samples"] = { **dict.fromkeys(VARIABLE_IDS, False), - context.variable: True, + variable: True, "formats": ["CSV"], "aggregation": previous_samples.get("aggregation", "point"), } @@ -245,9 +335,10 @@ def _describe(error: BaseException) -> str: def _write_record( context: EvaluationContext, + started_at: str, parameters: dict[str, float], nse: float | None, - station_nse: dict[str, float | None], + metrics: dict[str, StationMetrics], value: float, elapsed_seconds: float, error: str | None, @@ -258,13 +349,25 @@ def _write_record( 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_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, diff --git a/rubem/calibration/runner.py b/rubem/calibration/runner.py index 3d73fd21..8ffc5501 100644 --- a/rubem/calibration/runner.py +++ b/rubem/calibration/runner.py @@ -7,14 +7,20 @@ would accumulate the state of every run it made; ``spawn`` with ``max_tasks_per_child=1`` gives every evaluation a fresh interpreter. -The run directory receives three artifacts: ``evaluations.csv``, one row per -evaluation, ``result.json``, the best candidate and the settings that found it, -and ``-calibrated.json``, the calibrated configuration in the -format the input was written in. The rows of the table are ordered by the +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. It does not make the table itself reproducible: the ``id``, ``pid`` and -``elapsed_seconds`` columns differ between two runs of one seed, and from the +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`` @@ -22,10 +28,12 @@ raster in the user guide). Number of evaluations - The population of a ``sobol`` initialization is ``max(5, popsize * 8)`` - rounded up to the next power of two, since a Sobol' sequence is balanced - only over a power-of-two sample: with the default ``popsize`` of 15 that is - 128 members, so a calibration costs one model run per member and per + 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 @@ -42,9 +50,11 @@ 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 @@ -52,33 +62,44 @@ 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 -from .objective import Series, read_series +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, - bounds, - parameters_to_vector, - vector_to_parameters, - weights_constraint, + 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", @@ -87,6 +108,35 @@ ) """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.""" @@ -130,7 +180,9 @@ class CalibrationSettings: :type temp_dir: str | None :param init: Initial population: the name of a SciPy initialization - (``sobol`` by default) or an explicit ``(S, 8)`` array with ``S > 4``. + (``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. @@ -146,6 +198,24 @@ class CalibrationSettings: 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" @@ -160,6 +230,9 @@ class CalibrationSettings: 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) @@ -200,6 +273,14 @@ class CalibrationResult: :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] @@ -213,6 +294,8 @@ class CalibrationResult: evaluations_csv: Path result_json: Path calibrated_config: Path + best_series: Path + stations_csv: Path def _prepare_worker() -> None: @@ -257,7 +340,12 @@ class _Progress: generations: int = field(default=0) def __call__(self, intermediate_result) -> bool: - """Log the generation, the best objective so far and the records written. + """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 @@ -265,11 +353,15 @@ def __call__(self, intermediate_result) -> bool: ``False``. """ self.generations = int(getattr(intermediate_result, "nit", 0) or 0) - logger.info( - "Generation %d: best objective %.6g, %d evaluation(s) recorded.", + 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), - len(_record_paths(self.evaluations_dir)), + "n/a" if best_nse is None else f"{float(best_nse):.6f}", + len(records), ) return False @@ -301,8 +393,16 @@ def calibrate( :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, or if every - evaluation of the run failed. + 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. @@ -325,6 +425,18 @@ def calibrate( 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 @@ -333,19 +445,26 @@ def calibrate( 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=( - None - if aggregation is Aggregation.ZONES - else _station_ids(configuration.raster_files.sample_locations) - ), + station_ids=sampled_ids, ) - x0 = parameters_to_vector(configuration.calibration_parameters.model_dump()) + _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 @@ -357,6 +476,9 @@ def calibrate( "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 @@ -376,49 +498,70 @@ def calibrate( temp_dir=str(temp_dir), evaluations_dir=str(evaluations_dir), validate_input=False, + space=space, + stations=settings.stations, ) - planned_members = _population_size(settings) - logger.info( + 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).", - len(FREE_PARAMETERS), + space.dimension, planned_members, settings.maxiter, planned_members * (settings.maxiter + 1), workers, ) - - progress = _Progress(evaluations_dir) - executor = ProcessPoolExecutor( - max_workers=workers, - mp_context=multiprocessing.get_context("spawn"), - max_tasks_per_child=1, - initializer=_prepare_worker, - ) + _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: - optimum = differential_evolution( - functools.partial(evaluate, context=context), - bounds(), - strategy=settings.strategy, - maxiter=settings.maxiter, - popsize=settings.popsize, - mutation=settings.mutation, - recombination=settings.recombination, - rng=settings.seed, - polish=settings.polish, - init=settings.init, - x0=x0, - updating="deferred", - workers=executor.map, - constraints=weights_constraint(), - callback=progress, - ) - finally: - executor.shutdown(wait=True) + 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) - evaluations_csv = run_directory / EVALUATIONS_CSV _write_evaluations_csv(evaluations_csv, records) _report_failures(records) @@ -441,10 +584,26 @@ def calibrate( planned_members, ) - best_parameters = vector_to_parameters(optimum.x) + 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, @@ -452,17 +611,26 @@ def calibrate( 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), + }, ) - _write_calibrated_config(calibrated_config, configuration, file_v1, best_parameters) - logger.info( - "Calibration finished after %d evaluation(s) in %d generation(s): best objective %.6g.", + 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, @@ -476,6 +644,8 @@ def calibrate( evaluations_csv=evaluations_csv, result_json=result_json, calibrated_config=calibrated_config, + best_series=best_series, + stations_csv=stations_csv, ) @@ -494,6 +664,193 @@ def _require_scipy(): 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: @@ -634,16 +991,196 @@ def _check_observed_series( ) -def _population_size(settings: CalibrationSettings) -> int: +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 * 8)`` members, and ``sobol`` - rounds that up to the next power of two. + 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 * len(FREE_PARAMETERS)) + members = max(5, settings.popsize * dimension) if settings.init == "sobol": return 1 << (members - 1).bit_length() return members @@ -701,6 +1238,7 @@ def _write_evaluations_csv(path: Path, records: list[dict[str, Any]]) -> None: [ 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")), @@ -716,8 +1254,17 @@ def _cell(value) -> str: def _best_record(records: list[dict[str, Any]]) -> dict[str, Any] | None: - """Return the record of the lowest objective among those that were evaluated.""" - evaluated = [record for record in records if not record.get("error")] + """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"])) @@ -736,11 +1283,17 @@ def _write_result_json( 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.""" + """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, @@ -750,19 +1303,22 @@ def _write_result_json( "success": bool(optimum.success), "message": str(optimum.message), "population_size": population_size, - "settings": _settings_document(settings, workers, temp_dir), + "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, workers: int, temp_dir: Path + 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 and the temporary directory are the values the run - actually used, not the ``None`` the caller may have left them at, so that - the summary describes a calibration that can be repeated. + 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 { @@ -778,6 +1334,12 @@ def _settings_document( "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, } diff --git a/tests/unit/calibration/test_runner.py b/tests/unit/calibration/test_runner.py index 0d340619..856e65a5 100644 --- a/tests/unit/calibration/test_runner.py +++ b/tests/unit/calibration/test_runner.py @@ -10,6 +10,7 @@ import subprocess import sys from concurrent.futures import ProcessPoolExecutor +from concurrent.futures.process import BrokenProcessPool from pathlib import Path import numpy as np @@ -27,6 +28,8 @@ ) from rubem.calibration.runner import ( EVALUATION_COLUMNS, + OBSERVED_COLUMNS, + STATION_COLUMNS, CalibrationError, CalibrationSettings, _check_observed_series, @@ -110,6 +113,35 @@ def read_evaluations(path): 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) @@ -219,16 +251,110 @@ def test_the_temporary_directory_is_left_empty(self, calibrated): assert list(data.temp_dir.iterdir()) == [] @pytest.mark.unit - def test_the_run_directory_holds_the_three_artifacts(self, calibrated): + 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.""" @@ -328,23 +454,45 @@ def test_the_callback_takes_the_argument_scipy_passes_by_keyword(self): @pytest.mark.unit def test_the_callback_logs_the_generation_and_lets_the_search_go_on( - self, tmp_path, caplog, capfd + self, tmp_path, progress_log, capfd ): evaluations = tmp_path / "evaluations" evaluations.mkdir() (evaluations / "one.json").write_text("{}", encoding="utf8") - progress = _Progress(evaluations) + callback = _Progress(evaluations) - with caplog.at_level(logging.INFO, logger="rubem.calibration.runner"): - stop = progress(scipy_optimize.OptimizeResult(fun=12.5, nit=3)) + stop = callback(scipy_optimize.OptimizeResult(fun=12.5, nit=3)) assert stop is False, "the callback never halts the search" - assert progress.generations == 3 - assert "Generation 3" in caplog.text - assert "12.5" in caplog.text - assert "1 evaluation" in caplog.text + assert 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 @@ -529,6 +677,36 @@ def test_an_observed_station_the_run_does_not_sample_is_warned_about_and_ignored 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( @@ -650,15 +828,26 @@ def test_another_initialization_uses_the_members_it_asks_for(self): 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", "latinhypercube", "random"]) - def test_the_announced_population_is_the_one_scipy_builds(self, popsize, init): + @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. + # 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(), + bounds()[:dimension], init=init, popsize=popsize, maxiter=0, @@ -667,7 +856,7 @@ def test_the_announced_population_is_the_one_scipy_builds(self, popsize, init): ) assert len(optimum.population) == _population_size( - CalibrationSettings(popsize=popsize, init=init) + CalibrationSettings(popsize=popsize, init=init), dimension ) @pytest.mark.unit @@ -727,3 +916,318 @@ def test_importing_the_package_does_not_import_scipy(self, tmp_path): ) 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 index 301d1011..35948e6c 100644 --- a/tests/unit/calibration/test_worker.py +++ b/tests/unit/calibration/test_worker.py @@ -1,5 +1,6 @@ """One evaluation of a candidate: the model run, the record and the failure modes.""" +import datetime import json import os import shutil @@ -8,7 +9,12 @@ import pytest from rubem.api import Model -from rubem.calibration._worker import EvaluationContext, _derived_document, evaluate +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, @@ -39,7 +45,7 @@ def __init__(self, tmp_path): self.temp_dir.mkdir() self.evaluations_dir.mkdir() - def context(self, observed_path=None, variable="arn", spinup_steps=0): + 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, @@ -50,6 +56,7 @@ def context(self, observed_path=None, variable="arn", spinup_steps=0): temp_dir=str(self.temp_dir), evaluations_dir=str(self.evaluations_dir), validate_input=False, + stations=stations, ) def records(self): @@ -218,3 +225,159 @@ def test_a_variable_the_run_does_not_write_is_recorded_as_the_error(self, datase (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 fe9c9ba81bc997b038780df84333c9b7adf7dcea Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Thu, 17 Sep 2026 11:11:47 -0300 Subject: [PATCH 4/6] Expose the bounds, fixed parameters, stations and optimizer options on the command line --bound NAME=MIN:MAX and --fix NAME=VALUE, repeatable and refused with the expected spelling when malformed, --stations ID[,ID...], --init, --strategy and --polish/--no-polish reach the calibration settings; the closing summary lists every artifact of the run. --- rubem/cli.py | 160 +++++++++++++++++++++++++ tests/unit/test_cli_calibrate.py | 198 +++++++++++++++++++++++++++++++ 2 files changed, 358 insertions(+) diff --git a/rubem/cli.py b/rubem/cli.py index a1e0044e..eb402753 100644 --- a/rubem/cli.py +++ b/rubem/cli.py @@ -11,6 +11,7 @@ import logging import logging.config import logging.handlers +import math import sys import warnings from collections.abc import Sequence @@ -71,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__}") @@ -227,8 +344,39 @@ def calibrate( 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() @@ -241,11 +389,20 @@ def calibrate( "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 @@ -282,6 +439,9 @@ def calibrate( 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!") diff --git a/tests/unit/test_cli_calibrate.py b/tests/unit/test_cli_calibrate.py index 3fb9f4f0..5a3f8c57 100644 --- a/tests/unit/test_cli_calibrate.py +++ b/tests/unit/test_cli_calibrate.py @@ -54,6 +54,8 @@ def __init__(self, run_dir, best_nse=0.87): 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): @@ -97,6 +99,12 @@ def test_the_paths_and_the_documented_defaults( 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( @@ -206,6 +214,9 @@ def test_the_documented_progress_and_summary_are_printed( 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}" @@ -460,6 +471,13 @@ def test_the_help_documents_every_option(self, capsys, restore_logging): "--seed", "--workers", "--temp-dir", + "--bound", + "--fix", + "--stations", + "--init", + "--strategy", + "--polish", + "--no-polish", ): assert option in output, f"missing option: {option}\n{output}" @@ -477,3 +495,183 @@ def test_importing_the_command_line_does_not_need_scipy(self): 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 8590d136b7c60ec823cd62ce7a32364d0cba50ca Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Thu, 17 Sep 2026 11:11:47 -0300 Subject: [PATCH 5/6] Document the calibration options, artifacts and definitions --- doc/source/calibration.rst | 458 ++++++++++++++++++++++++++++++------- doc/source/userguide.rst | 79 +++++-- 2 files changed, 435 insertions(+), 102 deletions(-) diff --git a/doc/source/calibration.rst b/doc/source/calibration.rst index a8bb24cd..ded9d882 100644 --- a/doc/source/calibration.rst +++ b/doc/source/calibration.rst @@ -13,8 +13,8 @@ resulting time series at the stations and scores it against the observations. The options of the command are listed in the user guide, under :ref:`userguide:Running RUBEM`. This page documents what happens between them: the objective, the parameters and their bounds, the settings of the search and -how to budget a run, the inputs it needs, what it writes, and how it uses the -machine it runs on. +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 ------------------------- @@ -23,9 +23,10 @@ The model has nine calibration parameters, the ``CALIBRATION`` section of the configuration file: :math:`\alpha`, :math:`b`, :math:`w_1`, :math:`w_2`, :math:`w_3`, :math:`RCD`, :math:`f`, :math:`\alpha_{GW}` and :math:`x`. They are described in :ref:`overview:Calibration and Validation`. The search covers -eight of them; the slope factor weight :math:`w_3` is not searched but derived -from the other two weights, because the three weights of the potential runoff -coefficient must add up to 1. +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 @@ -71,17 +72,28 @@ Which values count `````````````````` A pair counts only when both of its members count. A value does not count when -it is not finite (``NaN`` and the infinities), when it is ``-9999``, the gap -marker of the model's own tables and of the usual observation tables, or when +it is not finite (``NaN`` and the infinities), when it is **negative**, or when it is at or above ``1e30``. -That last threshold is how the PCRaster missing value ``1e31`` is masked. The -value does not round-trip exactly through ``Float32`` — written and read back -it becomes ``9.999999...e30`` — so the mask rejects everything at or above +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 ```````````````````````` @@ -92,8 +104,21 @@ A station has no efficiency when fewer than two of its pairs count, or when all of its valid observations are equal: a constant record has no variability for the efficiency to explain and would divide by zero. Such a station is reported as having no efficiency and is left out of the average, instead of -dragging it. The efficiency of a candidate is the mean over the stations that -do have one. +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 @@ -102,7 +127,8 @@ common with the simulated steps that survive the spin-up, and one whose station ids match none of the ids of the sample locations raster (an observed station the configuration does not sample is reported once and ignored). With the ``zones`` aggregation the station ids are only known once a run has written -``zones_mapping.csv``, so that check is skipped. +: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 @@ -210,7 +236,8 @@ configuration the model accepts. 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. + 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: @@ -220,7 +247,7 @@ The ninth parameter is derived rather than searched: \[w_3 = 1 - \left(w_1 + w_2\right)\] -so the search is eight-dimensional under the linear constraint +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 @@ -231,11 +258,76 @@ 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. +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 ---------- @@ -243,22 +335,31 @@ The search Settings ```````` -The search is SciPy's ``differential_evolution`` with the following settings, -which the command does not expose because they describe the method rather than -the run: - -``strategy="best1exp"`` - The trial vector is built from the best member of the population plus a - scaled difference of two random members, with exponential crossover. A - best-based strategy exploits the best candidate found so far and converges - in fewer generations than the random-base ones [STORN1997]_, at some cost in - exploration, which the dithered mutation below offsets. Fewer generations is - what matters when one evaluation is a full simulation. - -``init="sobol"`` - The initial population is a Sobol' sequence over the bounds instead of a - random or Latin hypercube sample: it covers the eight-dimensional box far - more evenly, so the initial sample is not wasted on clusters and gaps. +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 @@ -271,11 +372,6 @@ the run: rather than of its parent. A high value suits parameters that interact, as the three weights and the runoff coefficients do. -``polish=False`` - No local refinement of the best candidate at the end. The local search would - spend model runs estimating gradients of an objective that is read back from - a simulated table, which does not provide them reliably. - ``updating="deferred"`` The whole generation is evaluated before the population is updated. This is what makes the generation a batch that can be spread over the worker @@ -289,6 +385,25 @@ 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 ````````````````````` @@ -298,10 +413,12 @@ The number of population members is not ``--popsize``: it is :label: calibration-members :nowrap: - \[members = \max\left(5,\ popsize \times 8\right)\] + \[members = \max\left(5,\ popsize \times N\right)\] -with 8 the number of free parameters, rounded up to the next power of two, -because a Sobol' sequence is balanced only over a power-of-two sample. Every +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 @@ -311,11 +428,16 @@ one more round of evaluations on top of ``--maxiter`` generations, so \[runs \leq members \times \left(maxiter + 1\right)\] -With the defaults, ``--popsize 15`` gives :math:`15 \times 8 = 120`, rounded up -to **128 members**, and ``--maxiter 100`` gives at most -:math:`128 \times 101 = 12\,928` model runs. A smaller search, ``--popsize 5``, -gives :math:`5 \times 8 = 40`, rounded up to **64 members**, so +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 @@ -323,7 +445,8 @@ 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. +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 @@ -331,12 +454,12 @@ 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). + 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. +evaluations and of generations it really spent. That number, and not the +announced one, is what ran. Inputs ------ @@ -345,7 +468,8 @@ The observed series ``````````````````` ``--observed`` takes one file, in either of two layouts, recognized from its -content. +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 @@ -359,10 +483,33 @@ in the first column and one value per station: 2;15.9;4.0;1.2 3;-9999;4.4;1.0 -The **PCRaster TSS layout** is accepted as well, both in its documented form (a -title line, the number of columns, one line per column name, then the rows) and -in the header-less form the model itself writes, whose stations are then -numbered ``1``, ``2``, ... in column order. +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:: @@ -375,6 +522,10 @@ numbered ``1``, ``2``, ... in column order. ``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 @@ -422,31 +573,73 @@ every variable. Artifacts of a run ------------------ -The three artifacts of the run directory, ``-o``/``--run-dir``: +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``, ``alpha``, - ``beta``, ``w_1``, ``w_2``, ``w_3``, ``rcd``, ``f``, ``alpha_gw``, ``x``, - ``nse``, ``objective``, ``elapsed_seconds`` and ``error``. ``id`` identifies - the evaluation and ``pid`` the process that made it; ``nse`` is empty for a - candidate that produced no efficiency; ``error`` is empty for a successful - evaluation, ``inadmissible`` for a candidate rejected without a run, and the - type and message of the exception for a run that failed. The rows are - ordered by the candidate they evaluated, not by the moment they were - written, so that the table does not depend on how the parallel workers - happened to finish. + 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``, and a ``settings`` object with - ``variable``, ``spinup_steps``, ``maxiter``, ``popsize``, ``seed``, - ``workers``, ``temp_dir``, ``init``, ``strategy``, ``mutation``, - ``recombination`` and ``polish``. The number of workers and the temporary - directory are the values the run resolved and used, not the defaults left - unset on the command line, so the summary describes a calibration that can - be repeated. + 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 @@ -457,9 +650,10 @@ The three artifacts of the run directory, ``-o``/``--run-dir``: Next to them, the :file:`evaluations/` directory holds one JSON record per evaluation, written by the worker that made it. A record carries what the CSV -row carries and, in addition, ``station_nse``, the efficiency of each -individual station — which is where to look when the mean hides one station -behaving differently from the others. +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:: @@ -468,6 +662,19 @@ behaving differently from the others. 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 -------------------------------- @@ -504,6 +711,91 @@ 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 --------------------------------------- @@ -524,7 +816,9 @@ to 1,344 runs, is well under an hour. 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. + 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 @@ -541,6 +835,12 @@ Shorten the simulation period first 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 @@ -549,7 +849,9 @@ Continue from a calibrated configuration with a different seed, a longer period or a larger population. Keep the run directory - :file:`evaluations.csv` and :file:`result.json` are the record of how a + :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 and the whole response surface the search sampled. They are - small, and they are what makes a published calibration reproducible. + 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/userguide.rst b/doc/source/userguide.rst index 222c7716..0bc20e9f 100644 --- a/doc/source/userguide.rst +++ b/doc/source/userguide.rst @@ -1182,38 +1182,66 @@ Use ``-h`` or ``--help`` to get a brief description of each command and its argu Calibrate the model parameters against an observed series. Options: - -c, --configfile Path to the configuration file (JSON). - [required] - --observed Observed series at the sample stations (CSV or - PCRaster TSS). [required] - -o, --run-dir Directory for the calibration artifacts. - [required] - --variable Output variable compared with the observed - series. [default: arn] - --spinup-steps Leading time steps excluded from the NSE. - [default: 0; x>=0] - --maxiter Maximum number of generations. [default: 100; - x>=1] - --popsize Population multiplier of the search. [default: - 15; x>=1] - --seed Seed of the differential evolution. - --workers Parallel evaluations (default: all cores but - one). [x>=1] - --temp-dir Parent of the per-evaluation output directories. - -h, --help Show this message and exit. + -c, --configfile Path to the configuration file (JSON). + [required] + --observed Observed series at the sample stations (CSV + or PCRaster TSS). [required] + -o, --run-dir Directory for the calibration artifacts. + [required] + --variable Output variable compared with the observed + series. [default: arn] + --spinup-steps Leading time steps excluded from the NSE. + [default: 0; x>=0] + --maxiter Maximum number of generations. [default: + 100; x>=1] + --popsize Population multiplier of the search. + [default: 15; x>=1] + --seed Seed of the differential evolution. + --workers Parallel evaluations (default: all cores but + one). [x>=1] + --temp-dir Parent of the per-evaluation output + directories. + --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 ``evaluations.csv``, ``result.json`` -and ``-calibrated.json`` into the ``--run-dir``. The method itself (the -objective, the parameter bounds, the derivation of ``w3``, how to budget the -number of evaluations and what a run needs from the configuration) is +"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 + $ 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 @@ -1222,6 +1250,9 @@ documented in :doc:`Calibration `. 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. From f0656eee1b13aa1a3c7c7e6aacf19d9ddeb6edff Mon Sep 17 00:00:00 2001 From: Gabriel Soares <70075435+soaressgabriel@users.noreply.github.com> Date: Thu, 17 Sep 2026 11:11:47 -0300 Subject: [PATCH 6/6] Record the calibration options and artifacts in the changelog --- doc/source/changelog.rst | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index da34c43b..a79279c8 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -13,6 +13,16 @@ The format follows `Keep a Changelog `__. Added ````` +- ``rubem calibrate`` gains per-run bounds (``--bound NAME=MIN:MAX``, narrowing the + settings ranges), fixed parameters that leave the search (``--fix NAME=VALUE``), + the stations that enter the objective (``--stations``), the initialization, + strategy and polishing options, per-station goodness-of-fit (pairs, means, + standard deviations, Pearson correlation, RMSE, NSE), ``observed.csv``, + ``stations.csv`` and ``best_.csv`` in the run directory, a + ``started_at`` column in the evaluations, progress lines on the terminal and + the consolidation of the evaluations of an interrupted search. Negative + observed values are gaps, and time series files without their header are + refused (`#348 `__). - Added the calibration page of the documentation (the objective, the free parameters and their bounds, the derivation of ``w3``, the differential evolution settings and the evaluation budget, the observed series formats,