diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1366150b..53b27abd 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -113,7 +113,7 @@ jobs: MATRIX_PYTHON: ${{ matrix.python }} shell: bash -el {0} - name: Install RUBEM with test dependencies - run: pip install -e ".[test,docs]" + run: pip install -e ".[test,docs,calibration]" shell: bash -el {0} - name: Run the test suite run: pytest -n auto -ra --cov=rubem --cov-report=xml --junitxml="TestResults-Rubem.xml" diff --git a/ci/dev-requirements.in b/ci/dev-requirements.in new file mode 100644 index 00000000..eec9fe2e --- /dev/null +++ b/ci/dev-requirements.in @@ -0,0 +1,18 @@ +# Direct development requirements audited by the pip-audit job. Refresh the +# lock with: +# uv pip compile ci/dev-requirements.in -o ci/dev-requirements.lock --python-version 3.12 +# An existing lock keeps its pins; add --upgrade to move every package forward. +build +codespell +pre-commit +pyfakefs +pytest-cov +pytest-mock +pytest-xdist +ruff +scipy>=1.16.2 +sphinx-autodoc-typehints +sphinx-rtd-theme +sphinxcontrib-email +sphinxcontrib-youtube +twine diff --git a/ci/dev-requirements.lock b/ci/dev-requirements.lock index 0e3042c7..00e304db 100644 --- a/ci/dev-requirements.lock +++ b/ci/dev-requirements.lock @@ -1,11 +1,11 @@ # This file was autogenerated by uv via the following command: -# uv pip compile dev-requirements.in -o dev-requirements.lock --python-version 3.12 +# uv pip compile ci/dev-requirements.in -o ci/dev-requirements.lock --python-version 3.12 alabaster==1.0.0 # via sphinx babel==2.18.0 # via sphinx build==1.5.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in certifi==2026.7.22 # via requests cffi==2.1.1 @@ -15,7 +15,7 @@ cfgv==3.5.0 charset-normalizer==3.5.1 # via requests codespell==2.4.3 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in coverage==7.15.4 # via pytest-cov cryptography==50.0.0 @@ -73,6 +73,8 @@ nh3==0.3.7 # via readme-renderer nodeenv==1.10.0 # via pre-commit +numpy==2.5.3 + # via scipy packaging==26.3 # via # build @@ -86,11 +88,11 @@ pluggy==1.6.0 # pytest # pytest-cov pre-commit==4.6.2 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in pycparser==3.0 # via cffi pyfakefs==6.2.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in pygments==2.21.0 # via # pytest @@ -101,16 +103,15 @@ pyproject-hooks==1.2.0 # via build pytest==9.1.1 # via - # -r dev-requirements.in # pytest-cov # pytest-mock # pytest-xdist pytest-cov==7.1.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in pytest-mock==3.15.1 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in pytest-xdist==3.8.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in python-discovery==1.5.2 # via virtualenv pyyaml==6.0.3 @@ -132,29 +133,30 @@ rich==15.0.0 roman-numerals==4.1.0 # via sphinx ruff==0.16.4 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in +scipy==1.18.1 + # via -r ci/dev-requirements.in secretstorage==3.5.0 # via keyring snowballstemmer==3.1.1 # via sphinx sphinx==9.1.0 # via - # -r dev-requirements.in # sphinx-autodoc-typehints # sphinx-rtd-theme # sphinxcontrib-email # sphinxcontrib-jquery # sphinxcontrib-youtube sphinx-autodoc-typehints==3.13.2 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in sphinx-rtd-theme==3.1.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in sphinxcontrib-applehelp==2.0.0 # via sphinx sphinxcontrib-devhelp==2.0.0 # via sphinx sphinxcontrib-email==0.3.6 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in sphinxcontrib-htmlhelp==2.1.0 # via sphinx sphinxcontrib-jquery==4.1 @@ -166,9 +168,9 @@ sphinxcontrib-qthelp==2.0.0 sphinxcontrib-serializinghtml==2.0.0 # via sphinx sphinxcontrib-youtube==1.5.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in twine==7.0.0 - # via -r dev-requirements.in + # via -r ci/dev-requirements.in urllib3==2.7.0 # via # id diff --git a/doc/source/api.rst b/doc/source/api.rst index 4ee8244f..b32d2c46 100644 --- a/doc/source/api.rst +++ b/doc/source/api.rst @@ -196,6 +196,37 @@ point guard above applies to it. A thread pool is not a substitute: threads that call ``Model.run`` share the state of PCRaster, with the consequences described above. +Calibration +----------- + +``rubem calibrate`` is built on this API: each of its workers calls +``Model.run()`` in its own process, on a configuration derived from the one +being calibrated. The same search is available from Python as +``rubem.calibration.runner.calibrate``, which takes the configuration file, the +observed series, the run directory and a ``CalibrationSettings``: + +.. code-block:: python + + from rubem.calibration.runner import CalibrationSettings, calibrate + + if __name__ == "__main__": + result = calibrate( + "config.json", + "observed.csv", + "calibration", + CalibrationSettings(seed=42, maxiter=20, popsize=5), + ) + print(result.best_parameters, result.best_nse) + +The entry point guard is required: the calibration starts its workers with the +``spawn`` method, which imports the main module of the caller in every one of +them. SciPy is needed as well, from the ``rubem[calibration]`` extra. + +``rubem.calibration`` is **not** part of the stable surface yet: only +``rubem.api`` is, and the calibration package may change without a minor +version bump until it is promoted. The method behind the command is documented +in :doc:`Calibration `. + Logging ------- diff --git a/doc/source/calibration.rst b/doc/source/calibration.rst new file mode 100644 index 00000000..e97978d2 --- /dev/null +++ b/doc/source/calibration.rst @@ -0,0 +1,887 @@ +Calibration +=========== + +.. role:: raw-html(raw) + :format: html + +The ``rubem calibrate`` command fits the calibration parameters of the model to +an observed streamflow series with a differential evolution search +[STORN1997]_. Every candidate the search proposes is a complete simulation: +the command builds the configuration of the candidate, runs it, samples the +resulting time series at the stations and scores it against the observations. + +The options of the command are listed in the user guide, under +:ref:`userguide:Running RUBEM`. This page documents what happens between them: +the objective, the parameters and their bounds, the settings of the search and +how to budget a run, the inputs it needs, what it writes, what it reports while +it runs, and how it uses the machine it runs on. + +What the calibration fits +------------------------- + +The model has nine calibration parameters, the ``CALIBRATION`` section of the +configuration file: :math:`\alpha`, :math:`b`, :math:`w_1`, :math:`w_2`, +:math:`w_3`, :math:`RCD`, :math:`f`, :math:`\alpha_{GW}` and :math:`x`. They +are described in :ref:`overview:Calibration and Validation`. The search covers +at most eight of them; the slope factor weight :math:`w_3` is not searched but +derived from the other two weights, because the three weights of the potential +runoff coefficient must add up to 1, and ``--fix`` takes further parameters out +of the search for one run. + +The observations are a streamflow series at the sample stations of the +configuration, one column per station. By default the simulated counterpart is +``arn``, the accumulated total runoff, in +:raw-html:`m3s-1`, which is the variable an observed +streamflow record is directly comparable with. ``--variable`` selects another +output variable, one of ``itp``, ``bfw``, ``srn``, ``eta``, ``lfw``, ``rec``, +``smc``, ``rnf`` and ``arn``; note that ``rnf`` is the total runoff of a single +cell, in millimeters, not a discharge, so an observed series compared with it +must be in the same quantity. + +Each evaluation runs the whole simulation period of the configuration. The +calibration never changes the period, the inputs or the aggregation: of what +the results depend on, it changes the nine parameters and nothing else. + +The objective function +---------------------- + +The efficiency +`````````````` + +The agreement between a simulated and an observed series at one station is the +Nash-Sutcliffe efficiency [MCCUEN2006]_: + +.. math:: + :label: calibration-nse + :nowrap: + + \[NSE = 1 - \frac{\sum_{t}{\left(Q_{sim,t} - Q_{obs,t}\right)^2}}{\sum_{t}{\left(Q_{obs,t} - \overline{Q_{obs}}\right)^2}}\] + +where: + +- :math:`Q_{sim,t}` – simulated value of the station at time step :math:`t`; +- :math:`Q_{obs,t}` – observed value of the station at time step :math:`t`; +- :math:`\overline{Q_{obs}}` – mean of the observed values used in the sum. + +Both sums run over the same pairs, and the mean is the mean of the observed +values of those pairs: a step dropped because the simulation has no value for +it is dropped from the mean as well, so numerator and denominator always +describe the same sample. + +Which values count +`````````````````` + +A pair counts only when both of its members count. A value does not count when +it is not finite (``NaN`` and the infinities), when it is **negative**, or when +it is at or above ``1e30``. + +Every quantity the two series carry is a flux or a storage and cannot be +negative, so a negative entry is not a measurement but a gap: the ``-9999`` the +model and the usual observation tables write, the ``-1`` some gauge records +use, or any other negative marker of their own. All of them are dropped alike, +and nothing has to be declared for it. Zero is a value and not a gap: a station +may well record no flow. + +The ``1e30`` threshold is how the PCRaster missing value ``1e31`` is masked. +The value does not round-trip exactly through ``Float32`` — written and read +back it becomes ``9.999999...e30`` — so the mask rejects everything at or above +``1e30`` instead of comparing with the constant. No streamflow or runoff of a +real basin comes anywhere near that magnitude. An empty cell, or a cell that is +not a number at all, is read as ``NaN`` and masked the same way. + +How many observed values the mask dropped is reported per station before the +search starts, on the terminal and in the ``dropped`` column of +:file:`observed.csv`, so that a gauge that is nearly empty over the simulated +window is seen before the model runs and not after them. + +Stations and the average +```````````````````````` + +The two series are aligned on the station ids they share and on the time steps +they share; a station or a step present on one side only is ignored. + +A station has no efficiency when fewer than two of its pairs count, or when all +of its valid observations are equal: a constant record has no variability +for the efficiency to explain and would divide by zero. Such a station is +reported as having no efficiency and is left out of the average, instead of +dragging it. + +The efficiency of a candidate is the mean over the **selected** stations that +do have one. Without ``--stations`` every shared station is selected. +``--stations 1,2,3`` names the ids the mean is taken over, and the stations it +leaves out are still sampled, still compared and still measured: they carry +``false`` in the ``in_selection`` column of :file:`observed.csv` and of +:file:`stations.csv`, and their goodness of fit is recorded by every +evaluation. That is the calibration and validation split — the search is +driven by some of the gauges, and the others say what the parameters it found +do at a gauge that had no say in them. + +A selection that names none of the comparable stations is refused before the +search starts, since every evaluation would fail on it; a selection that names +some unknown id keeps the ids that exist and reports the others once. + +A mismatch between the observations and the configuration is not a poor +candidate, and the two kinds of it that can be told before any model run are +refused before the search starts: an observed series with no time step in +common with the simulated steps that survive the spin-up, and one whose station +ids match none of the ids of the sample locations raster (an observed station +the configuration does not sample is reported once and ignored). With the +``zones`` aggregation the station ids are only known once a run has written +:file:`zones_mapping.csv`, so that check, and the check of the selected +stations, are made against the observed series alone. + +What cannot be told in advance surfaces inside the evaluations: a series whose +shared stations all lack an efficiency, because every observation is missing or +constant, fails each evaluation after its simulation has run. The evaluation is +then recorded with the message of the error in its ``error`` column and with +the objective ``1e30``, like any other failure, and the search carries on; the +calibration ends with an error only once no evaluation of the whole run has +succeeded. + +.. warning:: + + Such a failure is paid for in model runs before it is reported. The + per-generation log line is the early sign: a best objective of ``1e+30`` from + the first generation on means that every candidate failed. Check the observed + file against a run of the configuration first, over a few time steps, as + `Running on a cluster or a large machine`_ describes. + +The spin-up window +`````````````````` + +``--spinup-steps`` excludes the leading time steps of the simulation from the +comparison: every step whose number is at most ``--spinup-steps`` is dropped +from both series before they are aligned. The simulation still runs from its +first step — the point is to let the storages of the model fill before their +effect is scored, not to shorten the run. + +From the efficiency to the objective +```````````````````````````````````` + +The differential evolution minimizes, so the efficiency is turned into a cost: + +.. math:: + :label: calibration-objective + :nowrap: + + \[FO = 1000 \cdot \left(100 \cdot \left(1 - NSE\right)\right)^2\] + +A perfect simulation, :math:`NSE = 1`, gives :math:`FO = 0`; a simulation no +better than the mean of the observations, :math:`NSE = 0`, gives +:math:`FO = 10^7`; and the cost grows quadratically as the efficiency falls +further. + +A candidate that is never run — one the admissibility check rejects — and a +candidate whose run or evaluation fails are both recorded with the objective +``1e30``. The worst value a real run can produce is far below it (an efficiency +of :math:`-100`, already an absurd simulation, gives about +:math:`1.0 \times 10^{11}`), so a rejected or failed candidate always ranks +behind every candidate that was actually evaluated. The value is finite on +purpose: it sorts, it is written to JSON, and it appears in the table of +evaluations like any other objective. + +Free parameters and their bounds +-------------------------------- + +The eight searched parameters, in the order of the decision vector, with the +bounds the application settings declare (:file:`rubem/appsettings.json`, +section ``value_ranges.variables``). These are the same ranges the +configuration loader validates against, so any candidate inside the bounds is a +configuration the model accepts. + +.. list-table:: Searched parameters + :header-rows: 1 + :widths: 20 20 40 20 + + * - Configuration key + - Symbol + - Description + - Bounds + * - ``alpha`` + - :math:`\alpha` + - Interception parameter + - :math:`0.01 \leq \alpha \leq 10` + * - ``b`` + - :math:`b` + - Rainfall intensity coefficient + - :math:`0.01 \leq b \leq 1` + * - ``w_1`` + - :math:`w_1` + - Land use factor weight + - :math:`0 \leq w_1 \leq 1` + * - ``w_2`` + - :math:`w_2` + - Soil factor weight + - :math:`0 \leq w_2 \leq 1` + * - ``rcd`` + - :math:`RCD` + - Regional consecutive dryness level + - :math:`1 \leq RCD \leq 10` + * - ``f`` + - :math:`f` + - Flow direction factor + - :math:`0.01 \leq f \leq 1` + * - ``alpha_gw`` + - :math:`\alpha_{GW}` + - Baseflow recession coefficient + - :math:`0.01 \leq \alpha_{GW} \leq 1` + * - ``x`` + - :math:`x` + - Flow recession coefficient + - :math:`0 \leq x \leq 1` + +.. note:: + + The configuration file writes the rainfall intensity coefficient as ``b``, + and so does the calibrated configuration the run writes; + :file:`evaluations.csv` and :file:`result.json` spell it ``beta``, which is + the name the parameter carries inside the package. They are the same + parameter. ``--bound`` and ``--fix`` accept either spelling, as they accept + ``w1`` and ``w2`` for ``w_1`` and ``w_2``. + +The ninth parameter is derived rather than searched: + +.. math:: + :label: calibration-w3 + :nowrap: + + \[w_3 = 1 - \left(w_1 + w_2\right)\] + +so a search over both weights is under the linear constraint +:math:`w_1 + w_2 \leq 1`. The constraint is handed to the optimizer, which +checks it before the objective and never spends a model run on a candidate +whose slope factor weight would come out negative. The derived weight must also +lie inside its own range, :math:`0 \leq w_3 \leq 1`; a candidate for which it +does not is rejected without a run and recorded with the error +``inadmissible``. Both checks are guards rather than a part of the search: the +optimizer keeps its trial vectors inside the bounds and drops the ones that +violate the constraint, so it does not propose an inadmissible candidate of its +own. + +Narrowing a range for one run +````````````````````````````` + +``--bound NAME=MIN:MAX``, repeatable, searches one parameter in a range of the +caller's choosing instead of the range of the application settings: + +.. code-block:: console + + $ rubem calibrate -c project-config.json --observed observed.csv -o calibration \ + --bound alpha=1:10 --bound w_1=0.1:1 --bound rcd=1:5 + +A bound may only **narrow** the settings range, never widen it: its minimum +must lie at or above the minimum of the settings, its maximum at or below their +maximum, and the minimum below the maximum. Anything else is refused before the +search starts, because a candidate outside the settings range is a +configuration the loader would reject anyway. Bounding the derived ``w_3`` is +refused as well: bound ``w_1`` and ``w_2`` instead. So is bounding a parameter +that is fixed, which has no range left to search. + +Narrowing is what a known basin buys: the ranges of the settings are the ones +the model accepts at all, and a range that past work has already narrowed +concentrates the population where the answer is instead of spending members on +values that were never plausible. + +Keeping a parameter out of the search +````````````````````````````````````` + +``--fix NAME=VALUE``, repeatable, pins a parameter for the whole run: + +.. code-block:: console + + $ rubem calibrate -c project-config.json --observed observed.csv -o calibration --fix x=0 + +The coordinate **leaves the decision vector**: the search loses one dimension, +the optimizer never moves the parameter, and the population is no longer spent +exploring it. The pinned value is nonetheless one of the nine parameters every +candidate is run with, and it is in the best parameters, in the calibrated +configuration and under ``settings.fixed`` in :file:`result.json`. + +A fixed value must lie inside the range of the application settings, and at +least one parameter must be left free; both are checked before the search +starts. The fixed parameters, and any range a ``--bound`` narrowed, are logged +by the ``rubem.calibration`` logger when the run starts, so a logging +configuration that shows informational records replays exactly what was +searched. + +The two weights are the case where fixing changes more than the dimension: + +- **Both weights searched**: the pair is kept under the linear constraint + :math:`w_1 + w_2 \leq 1`. +- **One weight fixed**: there is no pair left to constrain, so the constraint + becomes a bound on the weight that is still free — its maximum drops to + :math:`1 - w_{fixed} - w_{3,min}`, which is what keeps the derived + :math:`w_3` inside its own range. If that leaves the free weight a single + value (``--fix w_1=1`` leaves :math:`w_2 = 0`), that weight is fixed at it + as well and leaves the decision vector, so the dimension and the budget + describe the search that runs; if it leaves no value at all, the run is + refused instead of searching an empty range. +- **Both weights fixed**: :math:`w_3` is a single number. It is derived once + and checked before the search; a pair that does not add up to a valid third + weight is refused. + +The starting point of the search, ``x0``, is the parameter set of the +configuration being calibrated. The optimizer places it in the initial +population, so the configuration's own parameters are always among the +candidates of the first generation and a calibration can never return something +worse than what it started from without the table showing why. A coordinate of +that starting point which a narrowed bound excludes is moved onto the bound it +crosses, and the move is logged: the run starts from the admissible point +closest to the configuration rather than ending before its first evaluation. + +The search +---------- + +Settings +```````` + +The search is SciPy's ``differential_evolution``. Three of the settings that +decide how it moves are on the command line, because a run may have to +reproduce the settings of an earlier calibration; the budget, ``--maxiter`` and +``--popsize``, is `The evaluation budget`_ and ``--seed`` is below; the rest +describe the method and are fixed. + +``--strategy`` (default ``best1exp``) + How a trial vector is built. With the default, it comes from the best member + of the population plus a scaled difference of two random members, with + exponential crossover. A best-based strategy exploits the best candidate + found so far and converges in fewer generations than the random-base ones + [STORN1997]_, at some cost in exploration, which the dithered mutation below + offsets. Fewer generations is what matters when one evaluation is a full + simulation. Every strategy name SciPy accepts is accepted here. + +``--init`` (default ``sobol``) + The initial population: ``sobol``, ``latinhypercube``, ``halton`` or + ``random``. The default is a Sobol' sequence over the bounds, which covers + the box far more evenly than a random sample, so the initial evaluations are + not wasted on clusters and gaps. It is also the one initialization that + resizes the population — see `The evaluation budget`_. + +``--polish`` / ``--no-polish`` (default ``--no-polish``) + Whether the best candidate is refined by a local search once the last + generation is over; see `Polishing`_. + +``mutation=(0.5, 1.0)`` + A pair, not a number: the mutation constant is drawn anew for every + generation between the two values (*dithering*). It widens the range of step + sizes the search takes without any tuning and reduces the chance of + stagnating in a local optimum. + +``recombination=0.7`` + The probability with which a trial vector keeps a component of the mutant + rather than of its parent. A high value suits parameters that interact, as + the three weights and the runoff coefficients do. + +``updating="deferred"`` + The whole generation is evaluated before the population is updated. This is + what makes the generation a batch that can be spread over the worker + processes; the immediate alternative is inherently serial. + +``--seed`` seeds the random generator of the search: the same seed gives the +same initial population and the same random draws of the mutation and of the +crossover. That alone does not make two calibrations comparable candidate by +candidate, because from the first generation on the selection reads the +objectives, so the candidates a run proposes depend on the values the +simulations returned. Two runs follow the same path only when the simulation +itself is reproducible — see `A fixed drainage network`_. + +Polishing +````````` + +Polishing is **off by default**, and on request with ``--polish``. When it is +on, the best candidate of the last generation is handed to a local search +(L-BFGS-B, or trust-constr when both weights are searched and the linear +constraint is therefore present), which keeps the result only if it improves +the objective. + +The reason not to pay for it by default is the price of a gradient here. The +local search has no derivatives to work from, so it estimates them by finite +differences: one model run per free parameter, plus the point itself, for every +step it takes, and those runs are made one at a time, outside the batch a +generation is spread over. On the full search one such step is nine serial +simulations, where a generation is 128 parallel ones. What it estimates is a gradient of an +objective read back from a simulated table, which the model does not provide +reliably at that scale. The evaluations polishing makes are recorded like any +other, so they appear in :file:`evaluations.csv` and in ``nfev``. + +The evaluation budget +````````````````````` + +The number of population members is not ``--popsize``: it is + +.. math:: + :label: calibration-members + :nowrap: + + \[members = \max\left(5,\ popsize \times N\right)\] + +with :math:`N` the number of free parameters — eight, minus one for every +``--fix``. With ``--init sobol``, the default, that count is then rounded up to +the next power of two, because a Sobol' sequence is balanced only over a +power-of-two sample; the other three initializations use it as it is. Every +generation evaluates one candidate per member, and the initial population is +one more round of evaluations on top of ``--maxiter`` generations, so + +.. math:: + :label: calibration-runs + :nowrap: + + \[runs \leq members \times \left(maxiter + 1\right)\] + +With the defaults and no fixed parameter, ``--popsize 15`` gives +:math:`15 \times 8 = 120`, rounded up to **128 members**, and ``--maxiter 100`` +gives at most :math:`128 \times 101 = 12\,928` model runs. The same run with +``--init latinhypercube`` keeps its **120 members** and at most +:math:`120 \times 101 = 12\,120` runs. A smaller search, ``--popsize 5``, gives +:math:`5 \times 8 = 40`, rounded up to **64 members**, so +``--popsize 5 --maxiter 20`` costs at most :math:`64 \times 21 = 1\,344` runs. +Fixing parameters shrinks :math:`N` and with it the population: four fixed +parameters leave :math:`N = 4`, so ``--popsize 15`` gives +:math:`15 \times 4 = 60`, rounded up to **64 members**. + +Two things make the real count lower. Members whose weights violate +:math:`w_1 + w_2 \leq 1` are skipped by the optimizer without a model run and +leave no record. And the search stops as soon as the objective values of the +population agree: SciPy compares their standard deviation with its default +tolerance of 1% of their mean, and skips the comparison while any member of +the population is still infeasible. A search that converges this way stops +before ``--maxiter`` is reached. Polishing works the other way and adds +evaluations after the last generation. + +The number of members and the resulting upper bound are logged when the +calibration starts, so the budget can be checked before the machine is +committed to it: + +.. code-block:: console + + Calibrating 8 free parameter(s) with 128 population member(s) per generation and at most 100 generation(s): up to 12928 model run(s), on 15 worker process(es). + +The size of the population the search actually built is recorded as +``population_size`` in :file:`result.json`, together with the number of +evaluations and of generations it really spent. That number, and not the +announced one, is what ran. + +Inputs +------ + +The observed series +``````````````````` + +``--observed`` takes one file, in either of two layouts, recognized from its +content: a first line carrying the ``;`` separator is read as the CSV table of +the model, anything else as a PCRaster time series. + +The **CSV layout** is the one the model writes for its own time series tables: +``;`` as the separator, a header whose first cell is ``0`` and whose remaining +cells are the station ids, and one row per time step carrying the step number +in the first column and one value per station: + +.. code-block:: console + + 0;1;2;3 + 1;12.4;3.1;0.8 + 2;15.9;4.0;1.2 + 3;-9999;4.4;1.0 + +The first cell of that header is the label of the step column: ``0`` in the +tables the model writes, or any non-numeric word in a table written by hand. A +first cell that is another number is a row of values, which means the file has +no header; it is refused, naming both layouts. + +The **PCRaster TSS layout** is its documented form: a title line such as +``timeseries scalar``, the number of columns, one line per column name +beginning with the time step column, and then one whitespace separated row per +time step: + +.. code-block:: console + + timeseries scalar + 4 + timestep + 1 + 2 + 3 + 1 12.4 3.1 0.8 + 2 15.9 4.0 1.2 + +The header is required. A time series file without it is refused rather than +read with columns numbered ``1``, ``2``, ...: the columns of a time series are +stations, nothing in a head-less file says which station a column belongs to, +and a numbering would quietly compare the observations of one gauge with the +simulation of another. Adding the header to a file that lacks one is a title, a +count and one line per column, and the refusal names the layout it expects. + +.. note:: + + The station ids must be the ids of the sample locations of the + configuration. The easiest way to get the layout and the ids right is to run + the configuration once and use the table the run wrote as the template of + the observed file: :file:`tss_arn.csv` for the default ``point`` + aggregation, :file:`tss_arn_subcatchment.csv` or :file:`tss_arn_zones.csv` + for the other two. With ``zones`` the ids are the renumbered columns + ``1..N``, and the correspondence with the original zone ids is written to + :file:`zones_mapping.csv` in the output directory. + +Gaps are written as they come: ``-9999``, any other negative marker, an empty +cell or a word that is not a number. The mask of `Which values count`_ is what +decides, so nothing has to be converted before a file is handed over. + +The configuration must name the raster the aggregation samples from — the +sample locations raster for ``point`` and ``subcatchment``, the zones raster +for ``zones``. A configuration that does not name the raster its aggregation +needs is refused before any run: there would be no station to compare. + +A fixed drainage network +```````````````````````` + +When the configuration does not provide a Local Drain Direction raster, the +model derives one from the digital elevation model with PCRaster's +``lddcreate``. On the real basin DEMs that operation is not deterministic: two +runs of the same configuration produce drainage networks that differ in a few +cells, and the accumulated runoff ``arn``, which is routed over the network, +therefore differs between them. The other eight output variables are computed +cell by cell and are unaffected. + +A calibration on ``arn`` without a fixed network would compare series routed +over different networks from one evaluation to the next, and the differences +between candidates would be partly noise. **Calibrating on** ``arn`` +**requires** ``RASTERS.ldd`` **to name a fixed LDD raster.** Generate it once, +with the DEM of the basin, and keep it with the other inputs: + +.. code-block:: python + + import pcraster as pcr + + pcr.setclone("/basins/ipojuca/input/maps/clone.map") + dem = pcr.readmap("/basins/ipojuca/input/maps/dem.map") + ldd = pcr.lddcreate(dem, 1e31, 1e31, 1e31, 1e31) + pcr.report(ldd, "/basins/ipojuca/input/maps/ldd.map") + +.. code-block:: json + + { + "RASTERS": { + "ldd": "/basins/ipojuca/input/maps/ldd.map", + }, + } + +The same applies beyond ``arn`` when the time series are aggregated over +subcatchments: the subcatchment of each sample location is delineated over the +drainage network, so the areas the series are averaged over move with it, for +every variable. + +Inputs the validation rejects +````````````````````````````` + +The inputs are validated **once**, when the calibration loads the +configuration; the worker processes never revalidate them, since they do not +change during the search. A blocking problem therefore stops the calibration +before the first model run, which is the point: a search is hundreds or +thousands of runs, and an input the validation refuses would spend all of them +on a result that cannot be trusted. + +``--allow-blocking-problems`` searches anyway. The checks still run and every +problem is still reported — the non-blocking ones as warnings, the blocking +ones as errors — and the search then starts instead of stopping. It is the +counterpart of ``rubem run --allow-blocking-problems``, and there is no ``-s`` +on ``calibrate`` for it to conflict with. + +The value is written to :file:`result.json`, under ``settings``, because it +changes what the numbers beside it mean: the parameters were fitted on inputs +the validation rejected, and nothing else in the artifacts would say so. + +.. warning:: + + The option does not repair anything. A raster the validation blocks on still + reaches the model at every one of the evaluations, and the efficiency the + search maximizes is measured on what that produced. Whether the defect + reaches the compared stations is worth answering before the machine is + committed: a single ``rubem run --allow-blocking-problems`` writes the time + series of the stations, and values that are finite and plausible there say + the search has something real to fit. + +Artifacts of a run +------------------ + +What the run directory, ``-o``/``--run-dir``, receives. The tables of the +stations use ``;`` as their separator, the one the model writes its own tables +with; :file:`evaluations.csv` is comma-separated. + +:file:`observed.csv` + Written **before** the search, one row per station of the observed series, + with the columns ``station``, ``in_selection``, ``pairs_in_window``, + ``dropped``, ``mean``, ``std``, ``min`` and ``max``. ``in_selection`` says + whether the station is one the objective averages over, which is every + station without ``--stations``; a station of the observed file the + configuration does not sample is listed here all the same, and never reaches + the objective. ``pairs_in_window`` is how many of the compared steps it + observes with a value the mask accepts, ``dropped`` how many values on those + steps the mask rejected as gaps, and the remaining columns summarize the + values that were kept. It is what a mean efficiency is + read against afterwards, and it costs nothing to look at first. + +:file:`evaluations.csv` + One row per evaluation, with the columns ``id``, ``pid``, ``started_at``, + ``alpha``, ``beta``, ``w_1``, ``w_2``, ``w_3``, ``rcd``, ``f``, + ``alpha_gw``, ``x``, ``nse``, ``objective``, ``elapsed_seconds`` and + ``error``. ``id`` identifies the evaluation, ``pid`` the process that made + it and ``started_at`` the wall-clock moment it began, in UTC; ``nse`` is + empty for a candidate that produced no efficiency; ``error`` is empty for a + successful evaluation, ``inadmissible`` for a candidate rejected without a + run, and the type and message of the exception for a run that failed. The + rows are ordered by the candidate they evaluated, not by the moment they + were written, so that the table does not depend on how the parallel workers + happened to finish — and ``started_at`` is what puts them back in the order + they were made in, which is the order a convergence curve is drawn in. + +:file:`stations.csv` + The goodness of fit of the best candidate, one row per station the two + series share, with the columns ``station``, ``in_selection``, ``pairs``, + ``mean_observed``, ``std_observed``, ``mean_simulated``, ``std_simulated``, + ``r`` (Pearson correlation), ``rmse`` and ``nse``. The statistics are the + ones the evaluation of that candidate recorded, not a recomputation, so the + table and the record always agree. A station outside the selection is + measured like any other and marked as such: this is where the validation + gauges are read. + +:file:`best_.csv` + The series behind the reported efficiency: the column ``step``, then + ``observed_`` and ``simulated_`` for every shared station, one row + per compared step. The steps and the stations are the ones the objective + compared, and a value either series does not offer is left as an empty cell. + The search keeps no output of its own — every evaluation writes into a + temporary directory that is removed again — so the winning candidate is run + once more, after the search, to produce this table. It is the last thing the + run does before the summary, so a failure of that final run is reported with + :file:`stations.csv` and the calibrated configuration already written and + :file:`result.json` not written at all; the search itself has finished and + its table is complete. + +:file:`result.json` + The summary: ``best_parameters`` (the nine parameters, the derived ``w_3`` + included), ``best_nse``, ``best_objective``, ``nfev`` and ``nit`` (the + evaluations and generations the search spent), ``success`` and ``message`` + from the optimizer, ``population_size``, an ``artifacts`` object naming + every file above, and a ``settings`` object with ``variable``, + ``spinup_steps``, ``maxiter``, ``popsize``, ``seed``, ``workers``, + ``temp_dir``, ``init``, ``strategy``, ``mutation``, ``recombination``, + ``polish``, ``bounds``, ``fixed`` and ``stations``. The workers, the + temporary directory and the bounds are the values the run resolved and used + — the bounds one entry per searched parameter, the settings ranges already + narrowed — and not the defaults left unset on the command line, so the + summary describes a calibration that can be repeated. + +:file:`-calibrated.json` + The calibrated configuration: the input document with the best parameters + in place, written in the format the input was written in, legacy or + format 1.0. Its paths are the ones the loader resolved, which are absolute + even when the input file wrote them relative to its own directory; moving + the file to another machine therefore means fixing the paths. + +Next to them, the :file:`evaluations/` directory holds one JSON record per +evaluation, written by the worker that made it. A record carries what the CSV +row carries and, in addition, ``station_metrics``, the statistics of +:file:`stations.csv` for every station of that evaluation, and ``station_nse``, +the efficiencies alone — which is where to look when the mean hides one station +behaving differently from the others, for a candidate that is not the best one. + +.. warning:: + + A run directory that already holds evaluation records of an earlier + calibration is refused before the search starts: consolidating two runs into + one table would mix them. Use an empty run directory for every calibration, + and keep the finished ones. + +An interrupted run +`````````````````` + +A search that is interrupted, or that dies, after its first evaluation still +leaves the evaluations it made: the records are already on disk, and they are +consolidated into :file:`evaluations.csv` on the way out, with a message saying +where they are. Together with the :file:`observed.csv` written before the +search and the records themselves, that is the whole history of the run. What +needs the finished search is not written: :file:`result.json`, +:file:`stations.csv`, :file:`best_.csv` and the calibrated +configuration. Since a run directory is never reused, an interrupted run is +read from its table and a new directory is used for the next attempt. + +Process model of the calibration +-------------------------------- + +The parent process loads the configuration once and validates its input files +once. The workers do not validate them again: they do not change during the +calibration, and re-reading every raster of the basin per evaluation would +dominate the cost. + +Every evaluation then runs in a **fresh interpreter**. PCRaster keeps its state +process-wide — ``setclone`` is global to the process and the raster memory of a +run is not reclaimed — so a worker reused for a second evaluation would carry +the state of the first. The pool is therefore started with the ``spawn`` method +and configured for one task per worker: a worker starts, evaluates one +candidate and exits. The generation is spread over as many such workers as +``--workers`` allows, which is also what makes the calibration parallel at all. + +Each worker builds the configuration of its candidate from the document of the +calibrated configuration: the nine parameters become the candidate's, the +output directory becomes a temporary directory of its own under ``--temp-dir``, +every raster series is disabled (an evaluation reads no raster of a previous +run, and writing them would dominate the cost of the run) and the time series +are reduced to the calibrated variable, as CSV. The aggregation the user +configured is kept, since it decides which areas the observed stations +correspond to. The temporary directory is removed when the evaluation ends, +whether it succeeded or not. + +The workers are silenced: a null handler on the root logger keeps the warnings +of a configuration that is not validated again out of the terminal, and +standard output is pointed at the null device so that the newline the PCRaster +framework writes when an interpreter exits is not printed once per evaluation. +Nothing is lost: a failure of an evaluation is recorded in that evaluation's +JSON record and in its row of the table, and the parent logs one summary at the +end saying how many evaluations failed and what the first failure was. One +failing candidate never ends a calibration; it is ranked behind every candidate +that was evaluated. + +What a run prints +````````````````` + +The progress of a calibration goes to the ``rubem.progress`` logger, which the +command line prints on standard output, exactly as it does for a simulation; a +calibration started from Python says nothing until its host configures that +logger. Four kinds of line are written, in this order: + +.. code-block:: console + + The calibration compares 216 time step(s); the observed series covers 216 of them, the step(s) 13 to 228. + Station 1: 198 observed value(s) on the compared steps, 18 dropped as gaps. + Station 2 (not in the objective): 204 observed value(s) on the compared steps, 12 dropped as gaps. + Calibrating 8 free parameter(s) with 128 population member(s) per generation and at most 100 generation(s): up to 12928 model run(s), on 15 worker process(es). + Generation 1: best objective 1.60757e+06, best NSE 0.601240, 256 evaluation(s) recorded. + Generation 2: best objective 4.21033e+05, best NSE 0.794900, 384 evaluation(s) recorded. + Calibration finished after 2048 evaluation(s) in 16 generation(s): best objective 160757, best NSE 0.873210. + +The coverage lines say how much of the compared window the observations offer, +per station, and which stations the objective leaves out. The budget line is +the arithmetic of `The evaluation budget`_. A generation line reports the best +objective the optimizer holds and, beside it, the best efficiency among the +records the workers have written so far, with how many evaluations have been +recorded; a best objective of ``1e+30`` means no candidate has been evaluated +successfully yet. The closing line repeats the search as SciPy reports it. + +Everything else — the fixed parameters, the narrowed ranges, the stations of +the objective, the warnings about a starting point moved onto a bound, and the +summary of the failed evaluations — goes to the ``rubem.calibration`` logger. +The command line shows its warnings; its informational records need a logging +configuration that asks for them. + +When a worker is killed +``````````````````````` + +A worker that disappears before it finishes its evaluation ends the +calibration, and the most common reason is the system killing it for lack of +memory: every worker holds the rasters of one model run. The message says so, +advises a lower ``--workers`` and names the table of the evaluations made up to +that point, which is written on the way out. + +Definitions worth knowing +------------------------- + +Five things about the method that are easier to read once than to infer from a +table of results. + +The weights are derived exactly + ``w_3`` is computed as :math:`1 - (w_1 + w_2)`, a single subtraction, so the + three weights of a candidate add up to 1 exactly and no candidate is ever + run with weights that do not. There is no tolerance and no penalty term: a + pair of weights that would leave ``w_3`` outside its range is not scored at + all, it is rejected before the model runs. + +A time step is a month + The model computes a monthly water balance, so the steps the efficiency is + computed on are the months of ``SIM_TIME``, and a window of 228 steps is + nineteen years. The steps carry the numbers the model gives them, counted + from the alignment month of the raster series rather than from the start of + the simulated period, so the first compared step need not be ``1``. + ``--spinup-steps`` drops steps by that number, and the coverage line the run + prints says which steps are left. + +A ``point`` station should mark one cell + The value a time series carries for a station is the average over all the + cells that carry that id in the sample locations raster. With the ``point`` + aggregation a station is meant to be a gauge, so its id should mark exactly + one cell; an id painted over several cells is compared as the mean of them, + which is a different quantity from the discharge of a gauge. The + ``subcatchment`` and ``zones`` aggregations are the ones that average an + area on purpose. + +The reported best was really run + The best candidate is one the search evaluated: its row is in + :file:`evaluations.csv`, with the objective the run produced. Nothing is + extrapolated, interpolated or reconstructed from the population, and + :file:`stations.csv` repeats the statistics of that very evaluation. + +The efficiency of the winner, station by station + :file:`stations.csv` has it, one row per station, with ``in_selection`` + separating the gauges that drove the search from the ones that only watched. + The same numbers are in the JSON record of that evaluation, under + ``station_metrics`` and ``station_nse``, and :file:`best_.csv` + holds the two series the numbers come from, ready to plot. + +Running on a cluster or a large machine +--------------------------------------- + +The cost of a calibration is the number of model runs times the cost of one +run, divided by the number of workers. Measure the second factor before +committing a machine: one ``rubem run`` of the configuration is an upper bound +on one evaluation, which writes no raster series and does not validate the +inputs again, and so costs a little less. On the Ipojuca basin, 228 monthly +time steps, an evaluation is of the order of half a minute. At 25 seconds +each, a default search on that basin, up to 12,928 runs on 15 workers, is on +the order of six hours; ``--popsize 5 --maxiter 20``, up +to 1,344 runs, is well under an hour. + +``--workers`` + Defaults to one less than the number of CPUs, which leaves a core for the + parent. On a shared machine, or under a batch scheduler, set it to the + number of cores the job was actually given rather than letting it read the + whole machine. Each worker holds one model run in memory, plus the parent, + so the memory of the job is roughly ``--workers`` times the memory of a + single ``rubem run`` of the same configuration; that, and not the core + count, is usually what limits the number of workers on a large basin, and a + job that asks for more workers than its memory allows is stopped by the + system killing one of them. + +``--temp-dir`` + The per-evaluation output directories are created here. Point it at fast + local storage of the node — the scratch directory of the job, not a shared + network filesystem: every evaluation creates a directory, writes a small + table into it and removes it again, and a network filesystem turns that into + the slowest part of the evaluation. The run directory itself may stay on + shared storage; it receives one small JSON record per evaluation. + +Shorten the simulation period first + An exploratory calibration over a few months of ``SIM_TIME``, with a small + ``--popsize`` and ``--maxiter``, costs a fraction of the full search and + still shows whether the observed series, the station ids and the spin-up are + set up as intended. Move to the full period only once a short run finishes + and produces a sensible efficiency. + +Fix and narrow what is already known + ``--fix`` removes a dimension from the search and ``--bound`` shrinks one. + Both make a given budget go further on the parameters that are still open, + and both are recorded in :file:`result.json`, so a run calibrated with part + of its parameters pinned still says which ones and at what value. + +Continue from a calibrated configuration + :file:`-calibrated.json` is an ordinary configuration file. Passing + it back to ``rubem calibrate`` as ``-c`` starts a second search from the + parameters the first one found, since the starting point of a search is the + parameter set of the configuration it is given — useful to refine a result + with a different seed, a longer period or a larger population. + +Keep the run directory + :file:`evaluations.csv`, :file:`observed.csv`, :file:`stations.csv`, + :file:`best_.csv` and :file:`result.json` are the record of how a + parameter set was obtained: the settings, the seed, the number of + evaluations, what the observations offered and the whole response surface + the search sampled. They are small, and they are what makes a published + calibration reproducible. diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index aa562d4b..58d2a07a 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -13,6 +13,24 @@ The format follows `Keep a Changelog `__. Added ````` +- Added ``rubem calibrate``: a differential evolution (SciPy) over the free + calibration parameters, with ``w3`` derived from the other two weights, + minimizing ``1000 (100 (1 - NSE))^2`` on the station series of one output + variable (``arn`` by default) against an observed series given in the table + layout of the model's time series or as a PCRaster time series with its + header. Negative, ``-9999`` and PCRaster missing values are gaps; a spin-up + window, per-run bounds (``--bound``), fixed parameters (``--fix``), the + stations of the objective (``--stations``) and the initialization, strategy + and polishing of the search are options. Every candidate runs in a fresh + spawned worker through the Python API, and the run directory receives + ``observed.csv``, ``evaluations.csv``, ``stations.csv``, + ``best_.csv``, ``result.json`` and ``-calibrated.json``, + with progress lines on the terminal. SciPy comes with the optional extra + ``rubem[calibration]``. The calibration page of the documentation describes + the method, the evaluation budget, the fixed drainage network a calibration + on ``arn`` needs and the process model; a ``dataset`` pytest marker with the + ``RUBEM_DATASET_DIR`` convention carries a smoke test on the Ipojuca basin + (`#343 `__). - Added ``rubem run --allow-blocking-problems``, which runs the model even when the input validation finds blocking problems: the same checks run, every problem is still reported (the non-blocking ones as warnings, the blocking ones as errors, followed by one error line @@ -22,7 +40,11 @@ Added not take it. ``ModelConfiguration(..., allow_blocking_problems=False)`` is the library counterpart and keeps the full list in ``ModelConfiguration.problems``; ``Model.from_file`` and ``Model.from_config`` of ``rubem.api`` take the same keyword and an isolated run - rebuilds the configuration with it + rebuilds the configuration with it. ``rubem calibrate --allow-blocking-problems``, and + ``CalibrationSettings(allow_blocking_problems=True)``, do the same for a calibration, + whose inputs are validated once before the search and never again by the workers; the + value is recorded under ``settings`` in :file:`result.json`, since the parameters it + reports were fitted on inputs the validation rejected (`#352 `__). - Added ``rubem.api``, the public Python surface: ``Model.from_file`` and ``Model.from_config`` load a configuration file or document, ``Model.run()`` diff --git a/doc/source/index.rst b/doc/source/index.rst index d2eff553..36615ffc 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -32,6 +32,7 @@ All development activity is coordinated via the `RUBEM Github page Path to the configuration file (JSON). + [required] + --observed Observed series at the sample stations (CSV + or PCRaster TSS). [required] + -o, --run-dir Directory for the calibration artifacts. + [required] + --variable Output variable compared with the observed + series. [default: arn] + --spinup-steps Leading time steps excluded from the NSE. + [default: 0; x>=0] + --maxiter Maximum number of generations. [default: + 100; x>=1] + --popsize Population multiplier of the search. + [default: 15; x>=1] + --seed Seed of the differential evolution. + --workers Parallel evaluations (default: all cores but + one). [x>=1] + --temp-dir Parent of the per-evaluation output + directories. + --bound Narrow the range of one parameter: + NAME=MIN:MAX. + --fix Keep one parameter out of the search: + NAME=VALUE. + --stations Station ids of the objective, comma- + separated (default: every shared station). + --init + Initial population of the search. [default: + sobol] + --strategy Differential evolution strategy. [default: + best1exp] + --polish / --no-polish Refine the best candidate with a local + search. [default: no-polish] + --allow-blocking-problems Search even if the input validation finds + blocking problems; they are logged as + errors instead of stopping the + calibration. + -h, --help Show this message and exit. + +``calibrate`` needs SciPy, which comes with ``pip install +"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. ``--allow-blocking-problems`` is the +counterpart of the option of ``run``: the inputs are validated once, before the +search, and the option lets a configuration whose inputs carry blocking +problems be searched on anyway. The method itself (the objective, the +parameter bounds, the derivation of ``w3``, how to budget the number of +evaluations, what the run prints and what it needs from the configuration) is +documented in :doc:`Calibration `. + +.. code-block:: console + + $ rubem calibrate -c project-config.json --observed observed.csv -o calibration \ + --spinup-steps 12 + Loading configuration and validating inputs... + Calibration started... + The calibration compares 216 time step(s); the observed series covers 216 of them, the step(s) 13 to 228. + Station 1: 198 observed value(s) on the compared steps, 18 dropped as gaps. + Calibrating 8 free parameter(s) with 128 population member(s) per generation and at most 100 generation(s): up to 12928 model run(s), on 15 worker process(es). + Generation 1: best objective 1.60757e+06, best NSE 0.601240, 256 evaluation(s) recorded. + Generation 16: best objective 160757, best NSE 0.873210, 2048 evaluation(s) recorded. + Calibration finished after 2048 evaluation(s) in 16 generation(s): best objective 160757, best NSE 0.873210. + Calibration finished successfully! + Best NSE: 0.873210 + Best objective: 160757 + Evaluations: 2048 + Generations: 16 + Evaluations table: /basins/ipojuca/calibration/evaluations.csv + Result: /basins/ipojuca/calibration/result.json + Calibrated configuration: /basins/ipojuca/calibration/project-config-calibrated.json + Observed summary: /basins/ipojuca/calibration/observed.csv + Stations table: /basins/ipojuca/calibration/stations.csv + Best candidate series: /basins/ipojuca/calibration/best_arn.csv + Use ``-V`` or ``--version`` to get the version of the RUBEM. .. code-block:: console diff --git a/pyproject.toml b/pyproject.toml index bc3a9341..051c5ef4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -66,8 +66,13 @@ preprocessing = [ "pykrige>=1.7", "scikit-gstat>=1.0", ] +# 1.16.2 is the first release with CPython 3.14 wheels for Linux and Windows, +# the range of interpreters this project supports. +calibration = [ + "scipy>=1.16.2", +] dev = [ - "rubem[test,docs,preprocessing]", + "rubem[test,docs,preprocessing,calibration]", "build", "codespell", "pre-commit", @@ -91,6 +96,7 @@ markers = [ "integration: mark a test as a integration test.", "exact: mark a test that reproduces the goldens byte-for-byte on the golden environment.", "paper: mark a conformity test derived from the published model formulation (journal supplement).", + "dataset: mark a test that needs the published datasets (RUBEM_DATASET_DIR).", ] [tool.coverage.run] diff --git a/rubem/_deps.py b/rubem/_deps.py index a27c0548..d6e846e1 100644 --- a/rubem/_deps.py +++ b/rubem/_deps.py @@ -69,3 +69,53 @@ def require_preprocessing_deps() -> None: "This preprocessing tool needs the optional dependencies " f"{', '.join(missing)}. Install them with 'pip install \"rubem[preprocessing]\"'." ) + + +_CALIBRATION_DEPENDENCIES = ("scipy",) + + +def missing_calibration_deps() -> list[str]: + """Return the optional calibration dependencies that are not importable. + + As in :func:`missing_runtime_deps`, a dependency whose lookup itself fails + counts as missing: the calibration would not get past it either. + + :return: The missing dependency names, in the order they are checked. + """ + missing = [] + for name in _CALIBRATION_DEPENDENCIES: + try: + found = importlib.util.find_spec(name) is not None + except (ImportError, ValueError): + found = False + if not found: + missing.append(name) + return missing + + +def calibration_deps_message(missing: list[str] | None = None) -> str: + """Return the installation guidance for the missing calibration dependencies. + + The command line raises it as ``SystemExit`` and the runner as + ``CalibrationError``, so both front ends give the same instructions. + + :param missing: The missing dependency names, see + :func:`missing_calibration_deps`. Defaults to ``None``, which names + every optional calibration dependency. + """ + names = list(_CALIBRATION_DEPENDENCIES) if missing is None else missing + plural = "dependencies" if len(names) > 1 else "dependency" + return ( + f"The calibration command needs the optional {plural} {', '.join(names)}. " + "Install it with 'pip install \"rubem[calibration]\"'." + ) + + +def require_calibration_deps() -> None: + """Fail fast with guidance when SciPy is not importable. + + :raises SystemExit: If any optional calibration dependency is missing. + """ + missing = missing_calibration_deps() + if missing: + raise SystemExit(calibration_deps_message(missing)) diff --git a/rubem/calibration/__init__.py b/rubem/calibration/__init__.py new file mode 100644 index 00000000..35a716c9 --- /dev/null +++ b/rubem/calibration/__init__.py @@ -0,0 +1,12 @@ +"""Automatic calibration of the model parameters against observed streamflow. + +The package holds the objective function (Nash-Sutcliffe efficiency over the +sample stations), the decision vector of the free calibration parameters, the +worker that evaluates one candidate, and the runner that drives SciPy's +differential evolution. + +Importing this package does not import SciPy: the optional dependency is +imported inside the functions that need it, so that reading a series or +computing an efficiency works in an installation without the +``rubem[calibration]`` extra. +""" diff --git a/rubem/calibration/_worker.py b/rubem/calibration/_worker.py new file mode 100644 index 00000000..b0d1e8b9 --- /dev/null +++ b/rubem/calibration/_worker.py @@ -0,0 +1,384 @@ +"""Evaluation of one candidate, the entry point the calibration workers run. + +The module lives at the top level of the package and keeps its entry point, +:func:`evaluate`, at module level, because the process pool of the calibration +starts its workers with the ``spawn`` method: the child imports this module by +name to find the function it has to call, so the function must be picklable by +reference and everything it needs must travel as plain data in an +:class:`EvaluationContext`. + +One evaluation builds the configuration of the candidate from the document of +the calibrated configuration, runs the model in the current process, reads the +time series it wrote and compares it with the observed series. It never calls +:meth:`rubem.api.Model.run_isolated`: the worker is already the isolation layer, +and an isolated run inside it would nest one process pool in another. + +:func:`evaluate` never raises: a candidate outside the admissible region is +rejected without a run, and any failure of the run or of the evaluation is +recorded and turned into +:data:`rubem.calibration.objective.INADMISSIBLE_OBJECTIVE`, so that one failing +candidate does not end the calibration. :func:`simulate_best`, which the search +calls once at the end to keep the series of the candidate it chose, does raise: +there is one named candidate and nothing to rank a failure behind. +""" + +import copy +import dataclasses +import json +import logging +import os +import shutil +import tempfile +import time +import uuid +from collections.abc import Sequence +from dataclasses import dataclass, field +from datetime import UTC, datetime +from pathlib import Path + +import numpy as np + +from ..configuration.model_configuration_file_v1 import VARIABLE_IDS +from .objective import ( + INADMISSIBLE_OBJECTIVE, + Series, + StationMetrics, + evaluate_series, + objective, + read_series, +) +from .parameters import DecisionSpace, decision_space + +logger = logging.getLogger(__name__) + +_INADMISSIBLE_ERROR = "inadmissible" +"""The error recorded for a candidate that was rejected without being run.""" + + +@dataclass(frozen=True) +class EvaluationContext: + """Everything an evaluation needs, as data that crosses to a worker process. + + :param document: The calibrated configuration as a format 1.0 document, + with its paths already anchored. + :type document: dict + + :param base_dir: Directory the relative paths of the document are anchored + on, ``None`` when they are absolute already. + :type base_dir: str | None + + :param variable: Id of the output variable the efficiency is computed on + (``arn``, the accumulated total runoff, by default). + :type variable: str + + :param observed: The observed series, read once by the parent. + :type observed: rubem.calibration.objective.Series + + :param spinup_steps: Number of initial time steps excluded from the + efficiency. + :type spinup_steps: int + + :param temp_dir: Directory the per-evaluation output directories are + created in and removed from. + :type temp_dir: str + + :param evaluations_dir: Directory the JSON record of every evaluation is + written to. + :type evaluations_dir: str + + :param validate_input: Whether every evaluation revalidates the input files + and their content. ``False`` by default: the parent validated them once + before the calibration started, and they do not change during it. + :type validate_input: bool + + :param space: The decision vector of the run, which turns a candidate into + the nine parameters of the model. Defaults to the whole search, without + fixed parameters and without narrowed bounds. + :type space: rubem.calibration.parameters.DecisionSpace + + :param stations: Ids of the stations the objective averages. ``None``, the + default, averages every station the two series share; the others are + measured either way. + :type stations: tuple[str, ...] | None + """ + + document: dict + base_dir: str | None + variable: str + observed: Series + spinup_steps: int + temp_dir: str + evaluations_dir: str + validate_input: bool = False + space: DecisionSpace = field(default_factory=decision_space) + stations: tuple[str, ...] | None = None + + +def evaluate(vector: Sequence[float] | np.ndarray, context: EvaluationContext) -> float: + """Run one candidate and return the objective value of its simulation. + + A candidate outside the admissible region is rejected before any model run + and recorded with the error ``inadmissible``. Otherwise the model is run in + the current process, into a temporary output directory that is removed + again, and the efficiency of the resulting time series is turned into the + objective. Every evaluation, successful or not, leaves one JSON record in + the evaluations directory of the context. + + :param vector: The candidate, the free values in the order of the + :attr:`~rubem.calibration.parameters.DecisionSpace.free_names` of the + decision space of the context. + :type vector: collections.abc.Sequence[float] | numpy.ndarray + + :param context: What the evaluation runs against. + :type context: EvaluationContext + + :return: The objective value, or + :data:`rubem.calibration.objective.INADMISSIBLE_OBJECTIVE` when the + candidate was rejected or its evaluation failed. + :rtype: float + """ + started_at = datetime.now(UTC).isoformat() + started = time.perf_counter() + try: + parameters = context.space.to_parameters(vector) + except ValueError as error: + _write_record( + context, started_at, {}, None, {}, INADMISSIBLE_OBJECTIVE, 0.0, _describe(error) + ) + return INADMISSIBLE_OBJECTIVE + + if not context.space.is_admissible(vector): + logger.debug("Rejecting an inadmissible candidate without running the model.") + _write_record( + context, + started_at, + parameters, + None, + {}, + INADMISSIBLE_OBJECTIVE, + time.perf_counter() - started, + _INADMISSIBLE_ERROR, + ) + return INADMISSIBLE_OBJECTIVE + + output_dir = tempfile.mkdtemp(prefix="evaluation-", dir=context.temp_dir) + try: + nse, metrics = _run_and_evaluate(context, parameters, output_dir) + except Exception as error: + # One failing candidate must not end the calibration: it is recorded and + # ranked behind every candidate that was actually evaluated. + logger.warning("The evaluation of a candidate failed: %s", error) + _write_record( + context, + started_at, + parameters, + None, + {}, + INADMISSIBLE_OBJECTIVE, + time.perf_counter() - started, + _describe(error), + ) + return INADMISSIBLE_OBJECTIVE + finally: + shutil.rmtree(output_dir, ignore_errors=True) + + value = objective(nse) + _write_record( + context, + started_at, + parameters, + nse, + metrics, + value, + time.perf_counter() - started, + None, + ) + return value + + +def simulate_best( + document: dict, + base_dir: str | None, + parameters: dict[str, float], + output_dir: str, + variable: str, +) -> str: + """Run one candidate and keep what it wrote. + + This is how the calibration obtains the series of the candidate it chose: + the run is the very run an evaluation makes, from the same derived + configuration, but the output directory is left in place so that the table + can be read afterwards. It lives at module level and takes plain data for + the same reason :func:`evaluate` does: the spawned workers import this + module by name to find the function they have to call. + + Unlike :func:`evaluate` it raises: the caller asked for one named + candidate's series and has nothing to rank a failure behind. + + :param document: The calibrated configuration as a format 1.0 document. + :type document: dict + + :param base_dir: Directory the relative paths of the document are anchored + on, ``None`` when they are absolute already. + :type base_dir: str | None + + :param parameters: The nine calibration parameters of the candidate. + :type parameters: dict[str, float] + + :param output_dir: Directory the run writes into. It is not removed. + :type output_dir: str + + :param variable: Id of the output variable the table is wanted for. + :type variable: str + + :return: The path of the CSV table the run wrote for ``variable``. + :rtype: str + + :raises RuntimeError: If the run wrote no CSV table for the variable. + """ + from ..api import Model + + derived = _document_for(document, variable, parameters, output_dir) + result = Model.from_config(derived, validate_input=False, base_dir=base_dir).run() + return str(_time_series_table(result, variable)) + + +def _run_and_evaluate( + context: EvaluationContext, + parameters: dict[str, float], + output_dir: str, +) -> tuple[float, dict[str, StationMetrics]]: + """Run the candidate and return its mean efficiency and the statistics per station.""" + from ..api import Model + + document = _derived_document(context, parameters, output_dir) + result = Model.from_config( + document, validate_input=context.validate_input, base_dir=context.base_dir + ).run() + simulated = read_series(_time_series_table(result, context.variable)) + return evaluate_series(simulated, context.observed, context.spinup_steps, context.stations) + + +def _time_series_table(result, variable: str) -> Path: + """Return the CSV table the run wrote for ``variable``. + + :raises RuntimeError: If the run reports no table for the variable. + """ + tables = result.time_series.get(variable, ()) + for path in tables: + if path.suffix.lower() == ".csv": + return path + raise RuntimeError( + f"The run wrote no CSV time series for the variable '{variable}'; it reports " + f"{', '.join(str(path) for path in tables) or 'no table'}." + ) + + +def _derived_document( + context: EvaluationContext, + parameters: dict[str, float], + output_dir: str, +) -> dict: + """Return the configuration document of one candidate of an evaluation.""" + return _document_for(context.document, context.variable, parameters, output_dir) + + +def _document_for( + base_document: dict, + variable: str, + parameters: dict[str, float], + output_dir: str, +) -> dict: + """Return the configuration document of one candidate. + + The calibrated document is copied and four things are changed: the + calibration parameters become the candidate's, the outputs go to the + directory of the run, every raster series is disabled (an evaluation reads + no raster, and writing them would dominate its cost) and the time series are + reduced to the calibrated variable, as CSV. The aggregation the user + configured is kept, since it decides which areas the observed stations + correspond to. + """ + document = copy.deepcopy(base_document) + document["model_calibration_parameters"] = { + "alpha": parameters["alpha"], + "b": parameters["beta"], + "w_1": parameters["w_1"], + "w_2": parameters["w_2"], + "w_3": parameters["w_3"], + "rcd": parameters["rcd"], + "f": parameters["f"], + "alpha_gw": parameters["alpha_gw"], + "x": parameters["x"], + } + output = document.setdefault("model_simulation_output", {}) + previous_rasters = output.get("raster_series", {}) + previous_samples = output.get("time_series_samples", {}) + output["dir_path"] = output_dir + output["raster_series"] = { + "no_data_value": previous_rasters.get("no_data_value", -9999), + "formats": [], + } + output["time_series_samples"] = { + **dict.fromkeys(VARIABLE_IDS, False), + variable: True, + "formats": ["CSV"], + "aggregation": previous_samples.get("aggregation", "point"), + } + return document + + +def _describe(error: BaseException) -> str: + """Return the recorded description of an exception: its type and its message.""" + return f"{type(error).__name__}: {error}" + + +def _write_record( + context: EvaluationContext, + started_at: str, + parameters: dict[str, float], + nse: float | None, + metrics: dict[str, StationMetrics], + value: float, + elapsed_seconds: float, + error: str | None, +) -> None: + """Write the JSON record of one evaluation, atomically. + + Every worker writes into the same directory while the parent may be reading + it, so a record is written to a temporary file next to its destination and + renamed onto it: a reader either does not see the record yet or sees it + whole. The name is a fresh UUID, so two workers never collide. + + The record carries the goodness of fit of every station under + ``station_metrics`` and, under ``station_nse``, the efficiencies alone, + taken from the very same statistics so that the two can never disagree. + ``started_at`` is the wall-clock moment the evaluation began, in UTC: the + records are consolidated in the order of the candidates they evaluated, and + this is what puts them back in the order they were made in. + """ + station_metrics = { + station: dataclasses.asdict(measured) for station, measured in metrics.items() + } + record = { + "id": uuid.uuid4().hex, + "pid": os.getpid(), + "started_at": started_at, + "parameters": parameters, + "nse": nse, + "station_nse": {station: measured["nse"] for station, measured in station_metrics.items()}, + "station_metrics": station_metrics, + "objective": value, + "elapsed_seconds": elapsed_seconds, + "error": error, + } + directory = Path(context.evaluations_dir) + target = directory / f"{record['id']}.json" + handle, temporary = tempfile.mkstemp(prefix=f".{target.name}.", suffix=".tmp", dir=directory) + try: + with os.fdopen(handle, "w", encoding="utf-8") as file: + json.dump(record, file) + Path(temporary).replace(target) + except BaseException: + Path(temporary).unlink(missing_ok=True) + raise diff --git a/rubem/calibration/objective.py b/rubem/calibration/objective.py new file mode 100644 index 00000000..2272fb49 --- /dev/null +++ b/rubem/calibration/objective.py @@ -0,0 +1,627 @@ +"""The calibration objective: Nash-Sutcliffe efficiency over the sample stations. + +The module reads the observed and the simulated time series, masks the values +that stand for a missing observation, computes the goodness of fit of each +station and turns the mean Nash-Sutcliffe efficiency (NSE) into the objective +the optimizer minimizes. + +Nothing here imports SciPy: the efficiency, the remaining statistics and the +objective are evaluated with NumPy alone. +""" + +import csv +import logging +from collections.abc import Sequence +from dataclasses import dataclass +from pathlib import Path + +import numpy as np + +from .._paths import PathInput, as_path + +logger = logging.getLogger(__name__) + +MISSING_VALUES = (-9999.0,) +"""Values written by the model and by the usual observation tables for a gap. + +``-9999`` is the one the model writes, and it is negative, so the sign rule of +:func:`valid_mask` would drop it on its own; it is named here because a table +that writes it is saying "gap" and not "an implausible discharge". +""" + +PCRASTER_MISSING = 1e31 +"""The PCRaster missing value. + +Written as ``Float32`` it does not round-trip exactly (``float(np.float32(1e31))`` +is ``9.999999...e30``), so the mask rejects everything at or above ``1e30`` +instead of comparing with this constant: no streamflow or runoff of a real +basin comes anywhere near that magnitude. +""" + +_MISSING_THRESHOLD = 1e30 + +INADMISSIBLE_OBJECTIVE = 1.0e30 +"""The objective of a candidate that is not run. + +The worst objective a run can produce is bounded by the efficiencies it can +reach: an NSE of -100, already an absurd simulation, gives +``1000 * (100 * 101) ** 2``, about ``1.0e11``. ``1.0e30`` is far above any such +value, so an inadmissible or failed candidate is always ranked behind every +candidate that was actually evaluated. It is finite so that it sorts, is +written to JSON and lands in the table of evaluations like any other objective. +""" + + +def valid_mask(values: np.ndarray) -> np.ndarray: + """Return the mask of the values that count as an observation. + + A value counts when it is finite, is not negative, is not one of + :data:`MISSING_VALUES` and is below ``1e30``, which covers the PCRaster + missing value in any of its round-trips. + + Every quantity the two series carry is a flux or a storage and cannot be + negative, so a negative entry is not a measurement: the observation tables + of the sample stations write the gaps of a gauge as ``-9999``, as ``-1`` or + as another negative marker of their own, and the mask drops all of them + alike. Zero is a value, not a gap: a station may well record no flow. + + :param values: Values of one station, in step order. + :type values: numpy.ndarray + + :return: A boolean array of the shape of ``values``. + :rtype: numpy.ndarray + """ + array = np.asarray(values, dtype=np.float64) + mask = np.isfinite(array) & (array < _MISSING_THRESHOLD) & (array >= 0.0) + for missing in MISSING_VALUES: + mask &= array != missing + return mask + + +def nash_sutcliffe(simulated: np.ndarray, observed: np.ndarray) -> float | None: + """Return the Nash-Sutcliffe efficiency of one station, or ``None``. + + The efficiency is computed over the pairs whose two members are valid, from + its definition:: + + NSE = 1 - sum((simulated - observed) ** 2) / sum((observed - mean(observed)) ** 2) + + Both sums run over the same pairs, and the mean is the mean of the observed + values of those pairs: a step dropped because the simulation has no value + for it is dropped from the mean as well, so that the numerator and the + denominator always describe the same sample. + + :param simulated: Simulated values, in step order. + :type simulated: numpy.ndarray + + :param observed: Observed values, aligned with ``simulated``. + :type observed: numpy.ndarray + + :return: The efficiency, or ``None`` when fewer than two pairs are valid or + when the valid observed values are all equal, which leaves the + denominator at zero and the efficiency undefined. + :rtype: float | None + + :raises ValueError: If the two series do not have the same length. + """ + simulated = np.asarray(simulated, dtype=np.float64) + observed = np.asarray(observed, dtype=np.float64) + if simulated.shape != observed.shape: + raise ValueError( + "The simulated and the observed series must have the same shape, got " + f"{simulated.shape} and {observed.shape}." + ) + + mask = valid_mask(simulated) & valid_mask(observed) + if int(np.count_nonzero(mask)) < 2: + return None + + simulated_valid = simulated[mask] + observed_valid = observed[mask] + # The zero-variance guard: a station whose observations are all equal has no + # variability for the efficiency to explain. Constancy is tested on the + # values themselves, since the sum of the squared deviations of a constant + # series is only approximately zero in floating point. + if np.all(observed_valid == observed_valid[0]): + return None + + residuals = float(np.sum((simulated_valid - observed_valid) ** 2)) + variability = float(np.sum((observed_valid - np.mean(observed_valid)) ** 2)) + return 1.0 - residuals / variability + + +@dataclass(frozen=True) +class StationMetrics: + """The goodness of fit of one station over the pairs both series have. + + Every field but :attr:`pairs` is ``None`` where the statistic is not + defined: a mean and a root mean squared error need one pair, a standard + deviation with ``ddof=1`` and an efficiency need two, and a correlation and + an efficiency need the series they divide by to vary. + + :param pairs: Number of time steps at which both series carry a value the + mask of :func:`valid_mask` accepts. + :type pairs: int + + :param mean_observed: Mean of the observed values of those pairs. + :type mean_observed: float | None + + :param std_observed: Standard deviation of the observed values of those + pairs, with one degree of freedom (``ddof=1``). + :type std_observed: float | None + + :param mean_simulated: Mean of the simulated values of those pairs. + :type mean_simulated: float | None + + :param std_simulated: Standard deviation of the simulated values of those + pairs, with one degree of freedom (``ddof=1``). + :type std_simulated: float | None + + :param r: Pearson correlation coefficient of the two series over those + pairs, ``None`` when either of them is constant. + :type r: float | None + + :param rmse: Root mean squared error over those pairs. + :type rmse: float | None + + :param nse: Nash-Sutcliffe efficiency, as :func:`nash_sutcliffe` defines it. + :type nse: float | None + """ + + pairs: int + mean_observed: float | None + std_observed: float | None + mean_simulated: float | None + std_simulated: float | None + r: float | None + rmse: float | None + nse: float | None + + +def station_metrics(simulated: np.ndarray, observed: np.ndarray) -> StationMetrics: + """Return the goodness of fit of one station. + + Every statistic is computed over the pairs whose two members pass + :func:`valid_mask`, the same sample :func:`nash_sutcliffe` uses, so the + efficiency and the statistics that explain it always describe one sample. + + :param simulated: Simulated values, in step order. + :type simulated: numpy.ndarray + + :param observed: Observed values, aligned with ``simulated``. + :type observed: numpy.ndarray + + :return: The statistics of the station. + :rtype: StationMetrics + + :raises ValueError: If the two series do not have the same length. + """ + simulated = np.asarray(simulated, dtype=np.float64) + observed = np.asarray(observed, dtype=np.float64) + if simulated.shape != observed.shape: + raise ValueError( + "The simulated and the observed series must have the same shape, got " + f"{simulated.shape} and {observed.shape}." + ) + + mask = valid_mask(simulated) & valid_mask(observed) + pairs = int(np.count_nonzero(mask)) + if pairs == 0: + return StationMetrics(0, None, None, None, None, None, None, None) + + simulated_valid = simulated[mask] + observed_valid = observed[mask] + mean_simulated = float(np.mean(simulated_valid)) + mean_observed = float(np.mean(observed_valid)) + rmse = float(np.sqrt(np.mean((simulated_valid - observed_valid) ** 2))) + if pairs < 2: + return StationMetrics(pairs, mean_observed, None, mean_simulated, None, None, rmse, None) + + return StationMetrics( + pairs=pairs, + mean_observed=mean_observed, + std_observed=float(np.std(observed_valid, ddof=1)), + mean_simulated=mean_simulated, + std_simulated=float(np.std(simulated_valid, ddof=1)), + r=_correlation(simulated_valid, observed_valid), + rmse=rmse, + nse=nash_sutcliffe(simulated_valid, observed_valid), + ) + + +def _correlation(simulated: np.ndarray, observed: np.ndarray) -> float | None: + """Return the Pearson correlation of two samples, or ``None``. + + The coefficient is undefined when either sample is constant, which leaves + its own spread at zero and the quotient without a denominator. Constancy is + tested on the values themselves, since the sum of the squared deviations of + a constant sample is only approximately zero in floating point. + """ + if np.all(simulated == simulated[0]) or np.all(observed == observed[0]): + return None + simulated_deviations = simulated - np.mean(simulated) + observed_deviations = observed - np.mean(observed) + spread = float(np.sqrt(np.sum(simulated_deviations**2) * np.sum(observed_deviations**2))) + return float(np.sum(simulated_deviations * observed_deviations)) / spread + + +def objective(nse: float) -> float: + """Return the objective value of an efficiency, ``1000 (100 (1 - NSE))^2``. + + The objective is zero for a perfect simulation and grows quadratically as + the efficiency falls, which is what the differential evolution minimizes. + + :param nse: The Nash-Sutcliffe efficiency of the candidate. + :type nse: float + + :return: The objective value. + :rtype: float + """ + return 1000.0 * (100.0 * (1.0 - nse)) ** 2 + + +@dataclass(frozen=True) +class Series: + """A time series table: the steps and the values of each station. + + :param steps: The time steps, one entry per row, as integers. + :type steps: numpy.ndarray + + :param stations: Station id to the values of that station, as ``float64`` + arrays aligned with :attr:`steps`. + :type stations: dict[str, numpy.ndarray] + """ + + steps: np.ndarray + stations: dict[str, np.ndarray] + + +def _as_value(token: str) -> float: + """Return ``token`` as a float, or ``NaN`` when it is empty or not a number. + + An entry a table leaves empty or writes in a form that is not a number is a + gap; it becomes ``NaN`` and the mask drops it, as it drops the gaps written + as ``-9999`` or as the PCRaster missing value. + """ + try: + return float(token) + except ValueError: + return float("nan") + + +def _as_step(token: str, path: Path, number: int) -> int: + """Return the time step written in the first column of a row. + + :raises ValueError: If the token is not a number, naming the file and the + row, since a table whose first column is a date or a label is a common + mistake and the bare conversion error does not say where it came from. + """ + try: + return int(float(token)) + except ValueError: + raise ValueError( + f"Row {number} of the time series file {path} does not start with a time " + f"step: {token!r} is not a number." + ) from None + + +def _build(ids: list[str], steps: list[int], rows: list[list[str]], path: Path) -> Series: + """Assemble a :class:`Series` from the parsed header, steps and value rows.""" + duplicated = sorted({station for station in ids if ids.count(station) > 1}) + if duplicated: + raise ValueError( + f"The time series file {path} has repeated station ids: {', '.join(duplicated)}." + ) + if not ids: + raise ValueError(f"The time series file {path} has no station column.") + if not rows: + raise ValueError(f"The time series file {path} has no data rows.") + + values = np.array( + [[_as_value(token) for token in row] for row in rows], + dtype=np.float64, + ) + return Series( + steps=np.asarray(steps, dtype=np.int64), + stations={station: values[:, column] for column, station in enumerate(ids)}, + ) + + +_ENCODING = "utf-8-sig" +"""The encoding both readers open a time series file with. + +It is UTF-8 with the byte order mark removed when the file carries one. The +mark is invisible and belongs to no cell, but it is written at the start of the +first line: left in place it becomes part of the first header cell, where it +turns the time step of a header-less table into a word and hides the refusal +that :func:`_check_csv_header` owes the reader. +""" + +_LAYOUTS = ( + "a time series file is either the table the model writes, whose first line is the " + "header '0;;...' and whose every further line is one time step, or a PCRaster " + "time series, whose header is a title line, the number of columns, and one line per " + "column name starting with the time step column" +) +"""The two layouts a time series file is read in, named by every refusal of a header.""" + + +def _read_csv(path: Path) -> Series: + """Read the CSV layout of the model's time series tables. + + The separator is ``;``, the header is ``0;;...`` and every row + carries the time step in the first column and one value per station. + + :raises ValueError: If the first record is a row of values instead of a + header, which would otherwise read the first time step as the station + ids and lose it from the series. + """ + with path.open(encoding=_ENCODING, newline="") as table: + records = [ + row for row in csv.reader(table, delimiter=";") if any(cell.strip() for cell in row) + ] + if not records: + raise ValueError(f"The time series file {path} is empty.") + + header = [cell.strip() for cell in records[0]] + _check_csv_header(header[0], path) + ids = header[1:] + steps: list[int] = [] + rows: list[list[str]] = [] + for number, record in enumerate(records[1:], 2): + cells = [cell.strip() for cell in record] + if len(cells) != len(header): + raise ValueError( + f"Row {number} of the time series file {path} has {len(cells)} column(s), " + f"the header gives {len(header)}." + ) + steps.append(_as_step(cells[0], path, number)) + rows.append(cells[1:]) + return _build(ids, steps, rows, path) + + +def _check_csv_header(label: str, path: Path) -> None: + """Refuse a ``;``-separated table whose first record is already a time step. + + The first cell of the header is the label of the step column, ``0`` in the + tables the model writes and any non-numeric word in a table written by + hand. A first cell that is another number is a row of values: the file has + no header, its first time step would be read as the station ids and the + remaining columns would be labelled with that step's values. + + :raises ValueError: If the first cell is a number other than ``0``. + """ + try: + step = float(label) + except ValueError: + return + if step == 0: + return + raise ValueError( + f"The time series file {path} has no header: its first line, {label!r}, is " + f"already a time step, and the station ids would be read from it. Add one: " + f"{_LAYOUTS}." + ) + + +def _read_tss(path: Path) -> Series: + """Read a PCRaster time series (``.tss``) file. + + The layout is a title line, a line with the number of columns, one line per + column name starting with the time step column, and then one whitespace + separated row per time step. + + :raises ValueError: If the header is not there. The columns of a time + series carry the ids of the stations they were sampled at, and nothing + in a header-less file says which station a column belongs to; numbering + the columns instead would label them with a guess. + """ + lines = [line for line in path.read_text(encoding=_ENCODING).splitlines() if line.strip()] + if not lines: + raise ValueError(f"The time series file {path} is empty.") + + columns = int(lines[1].strip()) if len(lines) > 1 and lines[1].strip().isdigit() else 0 + names = [line.strip() for line in lines[2 : 2 + columns]] + if columns < 2 or len(names) != columns: + raise ValueError( + f"The time series file {path} has no header: PCRaster time series files carry " + f"one, and its column names are the ids of the stations the columns were " + f"sampled at. Add it: {_LAYOUTS}." + ) + ids = names[1:] + data_start = 2 + columns + + rows: list[list[str]] = [] + steps: list[int] = [] + for number, line in enumerate(lines[data_start:], data_start + 1): + cells = line.split() + if len(cells) != len(ids) + 1: + raise ValueError( + f"Row {number} of the time series file {path} has {len(cells)} column(s), " + f"the column names give {len(ids) + 1}." + ) + steps.append(_as_step(cells[0], path, number)) + rows.append(cells[1:]) + return _build(ids, steps, rows, path) + + +def read_series(path: PathInput) -> Series: + """Read a time series table, in the model's CSV layout or as a PCRaster TSS. + + The layout is recognised from the content: a first line carrying the ``;`` + separator is the CSV table the model writes, anything else is read as a + PCRaster time series. Both layouts name their stations in a header, and a + file without one is refused rather than read with numbered columns: the + columns of a time series are stations, and a numbering would compare the + observations of one gauge with the simulation of another. A byte order mark + at the start of the file is dropped with the encoding, so that it cannot + become part of the first header cell. Missing entries are kept as they are + written; the mask of :func:`valid_mask` is what decides whether a value + counts. + + :param path: The file to read. + :type path: str | os.PathLike[str] + + :return: The parsed series. + :rtype: Series + + :raises FileNotFoundError: If the file does not exist. + :raises ValueError: If the file is empty, carries no header, has no station + column, repeats a station id or has a row that does not match its + header. + """ + file_path = as_path(path) + text = file_path.read_text(encoding=_ENCODING) + first_line = next((line for line in text.splitlines() if line.strip()), "") + if ";" in first_line: + return _read_csv(file_path) + return _read_tss(file_path) + + +@dataclass(frozen=True) +class Alignment: + """The steps and the stations two series are compared on. + + :param stations: The station ids both series carry, sorted. + :type stations: list[str] + + :param steps: The time steps both series carry after the spin-up, sorted. + :type steps: numpy.ndarray + + :param simulated_index: Rows of the simulated series, one per shared step. + :type simulated_index: numpy.ndarray + + :param observed_index: Rows of the observed series, one per shared step. + :type observed_index: numpy.ndarray + """ + + stations: list[str] + steps: np.ndarray + simulated_index: np.ndarray + observed_index: np.ndarray + + +def align_series(simulated: Series, observed: Series, spinup_steps: int = 0) -> Alignment: + """Return the steps and the stations two series have in common. + + The two series are aligned on the steps they share, after the steps up to + ``spinup_steps`` have been dropped from both, and on the station ids they + share: a station or a step present on one side only is ignored. This is the + one definition of the compared window, used both by the objective and by the + table of the best candidate, so that the two always describe one comparison. + + A step repeated in a table is aligned on its first row, since that is the + occurrence :func:`numpy.intersect1d` reports for a repeated value. + + :param simulated: The series the run wrote. + :type simulated: Series + + :param observed: The observed series. + :type observed: Series + + :param spinup_steps: Number of initial steps to exclude, defaults to ``0``. + :type spinup_steps: int, optional + + :return: The shared stations, the shared steps and the rows they are at. + :rtype: Alignment + + :raises ValueError: If the two series share no station id or no time step. + """ + shared_ids = sorted(set(simulated.stations) & set(observed.stations)) + if not shared_ids: + raise ValueError( + "The simulated and the observed series have no station in common: the " + f"simulation has {', '.join(sorted(simulated.stations)) or 'none'} and the " + f"observations have {', '.join(sorted(observed.stations)) or 'none'}." + ) + + simulated_steps = np.asarray(simulated.steps, dtype=np.int64) + observed_steps = np.asarray(observed.steps, dtype=np.int64) + simulated_kept = simulated_steps > spinup_steps + observed_kept = observed_steps > spinup_steps + shared_steps, simulated_rows, observed_rows = np.intersect1d( + simulated_steps[simulated_kept], + observed_steps[observed_kept], + return_indices=True, + ) + if shared_steps.size == 0: + raise ValueError( + "The simulated and the observed series have no time step in common after " + f"the {spinup_steps} spin-up step(s): the simulation has " + f"{simulated_steps[simulated_kept].tolist()} and the observations have " + f"{observed_steps[observed_kept].tolist()}." + ) + + return Alignment( + stations=shared_ids, + steps=shared_steps, + simulated_index=np.flatnonzero(simulated_kept)[simulated_rows], + observed_index=np.flatnonzero(observed_kept)[observed_rows], + ) + + +def evaluate_series( + simulated: Series, + observed: Series, + spinup_steps: int = 0, + stations: Sequence[str] | None = None, +) -> tuple[float, dict[str, StationMetrics]]: + """Return the mean efficiency of a run and the goodness of fit of each station. + + The two series are aligned by :func:`align_series`, and every station they + share is measured, whether or not it enters the objective. ``stations`` + names the ones the mean efficiency is taken over, which is how a run + calibrates on some of its gauges and reports the others as a validation; + without it every shared station enters the mean. The mean is taken over the + selected stations that have an efficiency, so that a station of zero + variance, reported as ``None``, does not drag the mean. + + :param simulated: The series the run wrote. + :type simulated: Series + + :param observed: The observed series. + :type observed: Series + + :param spinup_steps: Number of initial steps to exclude, defaults to ``0``. + :type spinup_steps: int, optional + + :param stations: Ids of the stations the mean is taken over. ``None``, the + default, takes it over every shared station. + :type stations: collections.abc.Sequence[str], optional + + :return: The mean efficiency of the selected stations and the statistics of + every shared station. + :rtype: tuple[float, dict[str, StationMetrics]] + + :raises ValueError: If the two series share no station id, share no time + step, if the selection names no shared station, or if no selected + station has an efficiency. + """ + alignment = align_series(simulated, observed, spinup_steps) + metrics = { + station: station_metrics( + simulated.stations[station][alignment.simulated_index], + observed.stations[station][alignment.observed_index], + ) + for station in alignment.stations + } + + if stations is None: + selected = list(alignment.stations) + else: + selected = [station for station in alignment.stations if station in set(stations)] + if not selected: + raise ValueError( + f"The station(s) {', '.join(stations) or 'none'} of the objective are " + "none of the stations the two series share " + f"({', '.join(alignment.stations)})." + ) + + values = [metrics[station].nse for station in selected if metrics[station].nse is not None] + if not values: + raise ValueError( + "No station has a Nash-Sutcliffe efficiency: every station of the objective " + f"({', '.join(selected)}) has fewer than two valid pairs or observations " + "without variance." + ) + return float(np.mean(values)), metrics diff --git a/rubem/calibration/parameters.py b/rubem/calibration/parameters.py new file mode 100644 index 00000000..78c16efd --- /dev/null +++ b/rubem/calibration/parameters.py @@ -0,0 +1,511 @@ +"""The decision vector of the calibration and its translation to model parameters. + +The differential evolution searches the free calibration parameters, ``alpha``, +``beta``, ``w_1``, ``w_2``, ``rcd``, ``f``, ``alpha_gw`` and ``x``; the slope +factor weight ``w_3`` is not searched but derived from the other two weights, +since the three weights of the weighted runoff coefficient must add up to 1. A +search over all eight is therefore eight-dimensional with the linear constraint +``w_1 + w_2 <= 1``, which keeps the derived weight inside its own range. + +:class:`DecisionSpace`, built by :func:`decision_space`, is that vector for one +run. It narrows the bounds of a parameter to a range of the user's choosing and +takes a parameter out of the vector altogether, at a value the user pins it to: +a fixed parameter keeps its value in the nine parameters every candidate is run +with, but the optimizer never moves it, and the dimension of the search falls by +one. Fixing one of the two searched weights replaces the linear constraint with +a narrower bound on the other, since only one of them is left to vary. + +The bounds a run starts from are the ones the application settings declare +(``appsettings.json``, section ``value_ranges.variables``), the same ranges +:class:`rubem.configuration.calibration_parameters.CalibrationParameters` +validates against, so a candidate the optimizer proposes inside the bounds is a +configuration the model accepts. An override may only narrow them, never widen +them, for the same reason. + +The module-level :func:`bounds`, :func:`vector_to_parameters`, +:func:`parameters_to_vector`, :func:`is_admissible` and +:func:`weights_constraint` are the default space, the one without fixed +parameters and without overridden bounds. + +Nothing here imports SciPy at module level; +:meth:`DecisionSpace.weights_constraint`, the only method that needs it, imports +it when it is called. +""" + +import logging +import math +from collections.abc import Mapping, Sequence +from dataclasses import dataclass +from typing import TYPE_CHECKING, cast + +import numpy as np + +from ..configuration._ranges import variable_range + +if TYPE_CHECKING: + from scipy.optimize import LinearConstraint + +logger = logging.getLogger(__name__) + +FREE_PARAMETERS = ("alpha", "beta", "w_1", "w_2", "rcd", "f", "alpha_gw", "x") +"""The parameters the optimizer searches, in the order of the decision vector.""" + +DERIVED_PARAMETER = "w_3" +"""The slope factor weight, derived from ``w_1`` and ``w_2`` instead of searched.""" + +CALIBRATION_PARAMETERS = ( + "alpha", + "beta", + "w_1", + "w_2", + "w_3", + "rcd", + "f", + "alpha_gw", + "x", +) +"""The nine calibration parameters of the model, in the order of the configuration.""" + +_ALIASES = {"b": "beta", "w1": "w_1", "w2": "w_2", "w3": "w_3"} +"""Spellings the legacy configuration file accepts for a calibration parameter.""" + + +@dataclass(frozen=True) +class DecisionSpace: + """The decision vector of one calibration: what is searched and between which bounds. + + The instances are plain data and are built by :func:`decision_space`, which + is what validates them; they travel to the worker processes inside the + evaluation context, so they hold nothing but names and numbers. + + :param free_names: The searched parameters, in the order of the decision + vector: the entries of :data:`FREE_PARAMETERS` that are not fixed. + :type free_names: tuple[str, ...] + + :param fixed: The parameters that are not searched, by name, with the value + every candidate carries. + :type fixed: dict[str, float] + + :param bounds: One ``(minimum, maximum)`` pair per free parameter, in the + order of :attr:`free_names`. + :type bounds: tuple[tuple[float, float], ...] + """ + + free_names: tuple[str, ...] + fixed: dict[str, float] + bounds: tuple[tuple[float, float], ...] + + @property + def dimension(self) -> int: + """Number of coordinates of the decision vector. + + :return: The number of free parameters. + :rtype: int + """ + return len(self.free_names) + + def to_parameters(self, vector: Sequence[float] | np.ndarray) -> dict[str, float]: + """Return the nine calibration parameters a decision vector stands for. + + The free values are read in the order of :attr:`free_names`, the fixed + ones are inserted at the value they were pinned to, and the slope factor + weight is derived as ``w_3 = 1 - w_1 - w_2``, so that the three weights + always add up to 1. + + :param vector: The free values, in the order of :attr:`free_names`. + :type vector: collections.abc.Sequence[float] | numpy.ndarray + + :return: The nine parameters, keyed by the names of + :class:`rubem.configuration.calibration_parameters.CalibrationParameters`. + :rtype: dict[str, float] + + :raises ValueError: If the vector does not have :attr:`dimension` entries. + """ + values = np.asarray(vector, dtype=np.float64).ravel() + if values.size != self.dimension: + raise ValueError( + f"The decision vector must have {self.dimension} entries " + f"({', '.join(self.free_names)}), got {values.size}." + ) + parameters = dict(self.fixed) + parameters.update( + {name: float(value) for name, value in zip(self.free_names, values, strict=True)} + ) + # The two weights are added first and the sum is subtracted once. Taking + # them away one at a time, ``1 - w_1 - w_2``, rounds twice: for 20 of the + # 101 pairs of two-decimal weights that add up to 1, among them + # ``w_1 = 0.33, w_2 = 0.67``, it lands one ulp below zero, and + # :class:`rubem.configuration.calibration_parameters.CalibrationParameters` + # then refuses the run. With a single subtraction the weight is exactly + # zero at the boundary and never negative while ``w_1 + w_2 <= 1``. + parameters[DERIVED_PARAMETER] = 1.0 - (parameters["w_1"] + parameters["w_2"]) + return {name: parameters[name] for name in CALIBRATION_PARAMETERS} + + def from_parameters(self, parameters: Mapping[str, float]) -> np.ndarray: + """Return the decision vector of a set of calibration parameters. + + This is how the configuration under calibration becomes the ``x0`` of the + differential evolution. Both spellings the configuration files use are + accepted, the canonical ``beta``, ``w_1``, ``w_2`` and the legacy ``b``, + ``w1``, ``w2``; the derived ``w_3`` and the fixed parameters are ignored, + since they are not searched. The straightforward way to obtain the + mapping is to dump the loaded + :class:`rubem.configuration.calibration_parameters.CalibrationParameters` + of a configuration, which has already normalised the spellings. + + :param parameters: The calibration parameters, by name. + :type parameters: collections.abc.Mapping[str, float] + + :return: The free values, in the order of :attr:`free_names`. + :rtype: numpy.ndarray + + :raises KeyError: If a free parameter is missing from the mapping. + """ + normalised = {_ALIASES.get(name, name): value for name, value in parameters.items()} + missing = [name for name in self.free_names if name not in normalised] + if missing: + raise KeyError( + f"The calibration parameters are missing {', '.join(missing)}; " + f"the decision vector needs {', '.join(self.free_names)}." + ) + return np.asarray([float(normalised[name]) for name in self.free_names], dtype=np.float64) + + def is_admissible(self, vector: Sequence[float] | np.ndarray) -> bool: + """Whether a decision vector describes a configuration the model accepts. + + A vector is admissible when every free value lies inside its bound, when + the two weights add up to at most 1 and when the derived ``w_3`` lies + inside its own bound. The derived weight is taken from + :meth:`to_parameters`, the very value the run would be given, so that the + guard can never disagree with the derivation and let through a candidate + :class:`rubem.configuration.calibration_parameters.CalibrationParameters` + would then refuse. + + :param vector: The free values, in the order of :attr:`free_names`. + :type vector: collections.abc.Sequence[float] | numpy.ndarray + + :return: ``True`` when the model can be run with this candidate. + :rtype: bool + """ + values = np.asarray(vector, dtype=np.float64).ravel() + if values.size != self.dimension: + return False + if not np.all(np.isfinite(values)): + return False + for value, (minimum, maximum) in zip(values, self.bounds, strict=True): + if not minimum <= value <= maximum: + return False + + parameters = self.to_parameters(values) + if parameters["w_1"] + parameters["w_2"] > 1.0: + return False + minimum, maximum = variable_range(DERIVED_PARAMETER) + return minimum <= parameters[DERIVED_PARAMETER] <= maximum + + def weights_constraint(self) -> "LinearConstraint | None": + """Return the linear constraint ``w_1 + w_2 <= 1`` on the decision vector. + + The differential evolution receives it so that it never spends a model + run on a candidate whose slope factor weight would come out negative: + SciPy evaluates the constraint before the objective and skips the members + that violate it. The constraint indexes the vector by position, so its + coefficients are placed at the positions the two weights occupy in + :attr:`free_names`. + + A search in which at most one of the two weights is free needs no + constraint: one fixed weight turns the sum into a bound on the other, + which :func:`decision_space` has already applied, and two fixed weights + leave nothing to constrain. + + :return: The constraint over the decision vector, or ``None`` when fewer + than two of the weights are searched. + :rtype: scipy.optimize.LinearConstraint | None + + :raises ImportError: If SciPy is not installed. + """ + positions = [ + self.free_names.index(name) for name in ("w_1", "w_2") if name in self.free_names + ] + if len(positions) < 2: + return None + + from scipy.optimize import LinearConstraint + + coefficients = np.zeros((1, self.dimension), dtype=np.float64) + for position in positions: + coefficients[0, position] = 1.0 + return LinearConstraint(coefficients, -np.inf, 1.0) + + +def decision_space( + fixed: Mapping[str, float] | None = None, + bounds: Mapping[str, tuple[float, float]] | None = None, +) -> DecisionSpace: + """Return the decision space of one calibration. + + Without arguments the space is the whole search: the eight free parameters + of :data:`FREE_PARAMETERS`, each between the bounds of the application + settings. + + ``fixed`` takes parameters out of the decision vector; their values stay in + the nine parameters every candidate is run with. ``bounds`` narrows the range + a free parameter is searched in. Both accept the legacy spellings ``b``, + ``w1`` and ``w2``. + + Fixing one of the two searched weights narrows the upper bound of the other + instead of constraining the pair, so that the derived ``w_3 = 1 - w_1 - w_2`` + cannot fall below its own minimum; when that narrowing leaves the other + weight a single value, it is fixed at that value as well and leaves the + decision vector. Fixing both of them derives ``w_3`` once and checks it here, + before the search starts. + + :param fixed: The parameters that are not searched, with the value every + candidate carries. Defaults to ``None``, no fixed parameter. + :type fixed: collections.abc.Mapping[str, float], optional + + :param bounds: The ``(minimum, maximum)`` range to search a parameter in, + by name. Defaults to ``None``, the ranges of the application settings. + :type bounds: collections.abc.Mapping[str, tuple[float, float]], optional + + :return: The decision vector of the run. + :rtype: DecisionSpace + + :raises ValueError: If a name is not a searched parameter or is the derived + ``w_3``, if a fixed value lies outside the range of the settings, if an + overridden bound is not a narrower range of the settings range, if a + fixed weight leaves the other one no admissible value, if two fixed + weights derive a ``w_3`` outside its range, or if every parameter is + fixed and nothing is left to search. + """ + fixed_values = { + _canonical(name, "fixed"): float(value) for name, value in (fixed or {}).items() + } + overrides = {_canonical(name, "bounded"): value for name, value in (bounds or {}).items()} + + settings_ranges = {name: variable_range(name) for name in FREE_PARAMETERS} + effective = dict(settings_ranges) + for name, override in overrides.items(): + if name in fixed_values: + raise ValueError( + f"'{name}' is fixed at {fixed_values[name]}, so a bound for it has nowhere " + "to apply: only the free parameters are searched." + ) + effective[name] = _narrowed(name, override, settings_ranges[name]) + + for name, value in fixed_values.items(): + minimum, maximum = settings_ranges[name] + if not math.isfinite(value) or not minimum <= value <= maximum: + raise ValueError( + f"'{name}' cannot be fixed at {value}: the value lies outside the range " + f"({minimum}, {maximum}) of the application settings." + ) + + free_names = tuple(name for name in FREE_PARAMETERS if name not in fixed_values) + if not free_names: + raise ValueError( + "Every calibration parameter is fixed, so the search has nothing to look for; " + "leave at least one of them free." + ) + + effective.update(_weight_bounds(fixed_values, effective)) + # A fixed weight can leave the other one a single admissible value (w_1 = 1 + # leaves w_2 = 0). A parameter with one value is not searched: it joins the + # fixed ones, so the dimension, the budget and the initial population describe + # the search that really runs, and no coordinate of zero width reaches SciPy. + for name in ("w_1", "w_2"): + if name in fixed_values: + continue + minimum, maximum = effective[name] + if minimum == maximum: + fixed_values[name] = minimum + free_names = tuple(free for free in free_names if free != name) + if not free_names: + raise ValueError( + "Every calibration parameter is fixed, so the search has nothing to look for; " + "leave at least one of them free." + ) + return DecisionSpace( + free_names=free_names, + fixed=fixed_values, + bounds=tuple(effective[name] for name in free_names), + ) + + +def _canonical(name: str, action: str) -> str: + """Return the canonical spelling of a parameter the caller named. + + :raises ValueError: If the name is the derived weight or is not a searched + parameter at all. + """ + canonical = _ALIASES.get(name, name) + if canonical == DERIVED_PARAMETER: + raise ValueError( + f"'{name}' cannot be {action}: the slope factor weight is derived from the " + "other two weights (w_3 = 1 - w_1 - w_2) and is not part of the decision vector; " + "fix or bound w_1 and w_2 instead." + ) + if canonical not in FREE_PARAMETERS: + raise ValueError( + f"'{name}' cannot be {action}: it is not a searched calibration parameter. " + f"The searched parameters are {', '.join(FREE_PARAMETERS)}." + ) + return canonical + + +def _narrowed( + name: str, override: Sequence[float], settings_range: tuple[float, float] +) -> tuple[float, float]: + """Return an overridden bound, checked against the range of the settings. + + :raises ValueError: If the override is not a range with a minimum below its + maximum, or if it is not contained in the range of the settings. + """ + try: + values = tuple(float(value) for value in override) + except TypeError: + # A single number instead of a pair: the same mistake as a pair of the + # wrong length, and it is answered with the same sentence. + raise ValueError( + f"The bound of '{name}' must be a (minimum, maximum) pair, got {override!r}." + ) from None + if len(values) != 2: + raise ValueError( + f"The bound of '{name}' must be a (minimum, maximum) pair, got {len(values)} value(s)." + ) + minimum, maximum = values + settings_minimum, settings_maximum = settings_range + if ( + not math.isfinite(minimum) + or not math.isfinite(maximum) + or minimum >= maximum + or minimum < settings_minimum + or maximum > settings_maximum + ): + raise ValueError( + f"The bound ({minimum}, {maximum}) of '{name}' is not a narrower range of its " + f"application settings range ({settings_minimum}, {settings_maximum}): a bound " + "narrows the settings range and its minimum lies below its maximum." + ) + return minimum, maximum + + +def _weight_bounds( + fixed_values: Mapping[str, float], effective: Mapping[str, tuple[float, float]] +) -> dict[str, tuple[float, float]]: + """Return the bounds the fixed weights impose on the weight that is left. + + With both weights searched the pair is handled by the linear constraint of + :meth:`DecisionSpace.weights_constraint` and nothing is narrowed here. With + one of them fixed the constraint would have a single term, which is a bound: + the free weight may not exceed what leaves the derived ``w_3`` at or above + its own minimum. With both of them fixed the derived weight is a number, and + it is checked here. + + :raises ValueError: If the free weight has no admissible value left, or if + two fixed weights derive a ``w_3`` outside its range. + """ + derived_minimum, derived_maximum = variable_range(DERIVED_PARAMETER) + weights = [name for name in ("w_1", "w_2") if name in fixed_values] + if not weights: + return {} + + if len(weights) == 2: + derived = 1.0 - (fixed_values["w_1"] + fixed_values["w_2"]) + if not derived_minimum <= derived <= derived_maximum: + raise ValueError( + f"Fixing w_1 at {fixed_values['w_1']} and w_2 at {fixed_values['w_2']} derives " + f"w_3 = {derived}, outside its range ({derived_minimum}, {derived_maximum}); " + "the three weights must add up to 1." + ) + return {} + + (pinned,) = weights + free_weight = "w_2" if pinned == "w_1" else "w_1" + minimum, maximum = effective[free_weight] + narrowed = min(maximum, 1.0 - fixed_values[pinned] - derived_minimum) + if narrowed < minimum: + raise ValueError( + f"Fixing {pinned} at {fixed_values[pinned]} leaves {free_weight} no value in its " + f"range ({minimum}, {maximum}): the derived w_3 = 1 - w_1 - w_2 would fall below " + f"{derived_minimum} for every one of them." + ) + return {free_weight: (minimum, narrowed)} + + +def bounds() -> list[tuple[float, float]]: + """Return the ``(minimum, maximum)`` bound of every free parameter. + + The bounds come from the application settings, in the order of + :data:`FREE_PARAMETERS`, which is the order of the decision vector of a + search without fixed parameters and without overridden bounds. + + :return: One ``(minimum, maximum)`` pair per free parameter. + :rtype: list[tuple[float, float]] + """ + return list(decision_space().bounds) + + +def vector_to_parameters(vector: Sequence[float] | np.ndarray) -> dict[str, float]: + """Return the nine calibration parameters a decision vector stands for. + + The vector is the one of the default decision space, the eight free values + in the order of :data:`FREE_PARAMETERS`; see + :meth:`DecisionSpace.to_parameters`. + + :param vector: The eight free values, in the order of :data:`FREE_PARAMETERS`. + :type vector: collections.abc.Sequence[float] | numpy.ndarray + + :return: The nine parameters, keyed by the names of + :class:`rubem.configuration.calibration_parameters.CalibrationParameters`. + :rtype: dict[str, float] + + :raises ValueError: If the vector does not have eight entries. + """ + return decision_space().to_parameters(vector) + + +def parameters_to_vector(parameters: Mapping[str, float]) -> np.ndarray: + """Return the decision vector of a set of calibration parameters. + + The vector is the one of the default decision space; see + :meth:`DecisionSpace.from_parameters`. + + :param parameters: The calibration parameters, by name. + :type parameters: collections.abc.Mapping[str, float] + + :return: The eight free values, in the order of :data:`FREE_PARAMETERS`. + :rtype: numpy.ndarray + + :raises KeyError: If a free parameter is missing from the mapping. + """ + return decision_space().from_parameters(parameters) + + +def is_admissible(vector: Sequence[float] | np.ndarray) -> bool: + """Whether a decision vector describes a configuration the model accepts. + + The vector is the one of the default decision space; see + :meth:`DecisionSpace.is_admissible`. + + :param vector: The eight free values, in the order of :data:`FREE_PARAMETERS`. + :type vector: collections.abc.Sequence[float] | numpy.ndarray + + :return: ``True`` when the model can be run with this candidate. + :rtype: bool + """ + return decision_space().is_admissible(vector) + + +def weights_constraint() -> "LinearConstraint": + """Return the linear constraint ``w_1 + w_2 <= 1`` on the decision vector. + + The constraint is the one of the default decision space, in which both + weights are searched; see :meth:`DecisionSpace.weights_constraint`. + + :return: The constraint, over the eight entries of the decision vector. + :rtype: scipy.optimize.LinearConstraint + + :raises ImportError: If SciPy is not installed. + """ + # Both weights are searched in the default space, so the constraint is there. + return cast("LinearConstraint", decision_space().weights_constraint()) diff --git a/rubem/calibration/runner.py b/rubem/calibration/runner.py new file mode 100644 index 00000000..e7610e5e --- /dev/null +++ b/rubem/calibration/runner.py @@ -0,0 +1,1397 @@ +"""The calibration itself: a differential evolution over the free parameters. + +:func:`calibrate` loads and validates the configuration once, reads the observed +series, and hands SciPy's ``differential_evolution`` a process pool whose workers +run one model evaluation each and exit. The pool is the isolation layer: PCRaster +keeps its clone and its raster memory process-wide, so a worker that is reused +would accumulate the state of every run it made; ``spawn`` with +``max_tasks_per_child=1`` gives every evaluation a fresh interpreter. + +The run directory receives, before the search, ``observed.csv``, the coverage +and the summary statistics of the observations over the compared steps, and +after it ``evaluations.csv``, one row per evaluation, ``stations.csv``, the +goodness of fit of the best candidate at every station, ``best_.csv``, +the observed and the simulated series of that candidate side by side, +``result.json``, the best candidate and the settings that found it, and +``-calibrated.json``, the calibrated configuration in the format +the input was written in. The rows of the table of evaluations are ordered by the +candidate they evaluated and not by the moment they were written, so that the +order of the rows does not depend on the order in which the workers happened to +finish; the ``started_at`` column is what puts them back in the order they were +made in. It does not make the table itself reproducible: the ``id``, ``pid``, +``started_at`` and ``elapsed_seconds`` columns differ between two runs of one +seed, and from the +first generation on the selection reads the objective values, so two runs +evaluate the same candidates only when the simulation itself is reproducible +(on the real basins that requires a fixed ``RASTERS.ldd``, since ``lddcreate`` +does not always derive the same directions twice; see the note on the LDD +raster in the user guide). + +Number of evaluations + The population of a ``sobol`` initialization is ``max(5, popsize * N)``, + with ``N`` the number of free parameters, rounded up to the next power of + two, since a Sobol' sequence is balanced only over a power-of-two sample: + with the default ``popsize`` of 15 and the eight free parameters of a search + without fixed ones that is 128 members, so a calibration costs one model run + per member and per + generation, ``128 * (maxiter + 1)`` runs at most. The number is logged when + the calibration starts so that the budget can be checked before it runs. + Members whose weights violate ``w_1 + w_2 <= 1`` are skipped by SciPy + without a model run, so the real count is usually lower. + +SciPy is an optional dependency, imported inside :func:`calibrate`; importing +this module without it works and raises nothing. +""" + +import csv +import functools +import json +import logging +import math +import multiprocessing +import os +import shutil +import sys +import tempfile +from concurrent.futures import ProcessPoolExecutor +from concurrent.futures.process import BrokenProcessPool +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any + +import numpy as np + +from .._paths import PathInput, as_path +from ..configuration._ranges import variable_range +from ..configuration.model_configuration_file_v1 import ( + VARIABLE_IDS, + Aggregation, + ModelConfigurationFileV1, +) +from ._worker import _INADMISSIBLE_ERROR, EvaluationContext, evaluate, simulate_best +from .objective import Series, align_series, read_series, valid_mask +from .parameters import ( + CALIBRATION_PARAMETERS, + FREE_PARAMETERS, + DecisionSpace, + decision_space, +) + +logger = logging.getLogger(__name__) + +progress = logging.getLogger("rubem.progress") +"""The logger of what a person watching a calibration wants to see. + +The budget, one line per generation and the closing summary go here; the module +logger keeps the diagnostics, the warnings and the refusals. The command line +routes this logger to standard output, as it does for a simulation, so that an +embedded calibration stays silent unless its host asks for the progress. +""" + +EVALUATIONS_DIRNAME = "evaluations" +"""Name of the directory under the run directory that collects the JSON records.""" + +EVALUATIONS_CSV = "evaluations.csv" +RESULT_JSON = "result.json" +OBSERVED_CSV = "observed.csv" +STATIONS_CSV = "stations.csv" + +EVALUATION_COLUMNS = ( + "id", + "pid", + "started_at", + *CALIBRATION_PARAMETERS, + "nse", + "objective", + "elapsed_seconds", + "error", +) +"""The columns of ``evaluations.csv``, one row per evaluation.""" + +OBSERVED_COLUMNS = ( + "station", + "in_selection", + "pairs_in_window", + "dropped", + "mean", + "std", + "min", + "max", +) +"""The columns of ``observed.csv``, one row per station of the observed series.""" + +STATION_COLUMNS = ( + "station", + "in_selection", + "pairs", + "mean_observed", + "std_observed", + "mean_simulated", + "std_simulated", + "r", + "rmse", + "nse", +) +"""The columns of ``stations.csv``, one row per station the best candidate compared.""" + +STATION_SEPARATOR = ";" +"""The separator of the per-station tables, the one the model writes its own tables with.""" + + +class CalibrationError(RuntimeError): + """A calibration that could not start, or that produced no usable evaluation.""" + + +@dataclass(frozen=True) +class CalibrationSettings: + """How the differential evolution is run. + + The differential evolution defaults follow the method the calibration of the + model was designed around; the remaining fields are the knobs the command + line exposes. + + :param variable: Id of the output variable compared with the observations. + ``arn``, the accumulated total runoff in m3/s, by default; ``rnf`` is + the total runoff per cell in mm. + :type variable: str + + :param spinup_steps: Number of initial time steps excluded from the + efficiency. + :type spinup_steps: int + + :param maxiter: Maximum number of generations. + :type maxiter: int + + :param popsize: Population multiplier; see the module docstring for the + number of members it produces. + :type popsize: int + + :param seed: Seed of the random generator of the search, for a reproducible + calibration. ``None`` leaves the generator unseeded. + :type seed: int | None + + :param workers: Number of worker processes. ``None``, the default, resolves + at call time to one less than the number of CPUs, at least one, which + leaves a core for the parent. + :type workers: int | None + + :param temp_dir: Directory the per-evaluation output directories are + created in. ``None`` uses the temporary directory of the system. + :type temp_dir: str | None + + :param init: Initial population: the name of a SciPy initialization + (``sobol`` by default) or an explicit ``(S, N)`` array with ``S > 4`` + members and one column per free parameter, which is one column fewer + for every parameter :attr:`fixed` takes out of the decision vector. + :type init: str | numpy.ndarray + + :param strategy: Differential evolution strategy. + :type strategy: str + + :param mutation: Mutation constant, a ``(min, max)`` pair for dithering. + :type mutation: tuple[float, float] + + :param recombination: Crossover probability. + :type recombination: float + + :param polish: Whether a local search refines the best candidate. ``False``: + the local search would spend model runs on gradients the objective, a + simulation read from a table, does not provide reliably. + :type polish: bool + + :param bounds: The ``(minimum, maximum)`` range to search a parameter in, by + name; it narrows the range of the application settings, which is what a + parameter left out of the mapping is searched in. ``None``, the default, + searches every parameter in the range of the settings. + :type bounds: dict[str, tuple[float, float]] | None + + :param fixed: The parameters that are not searched, by name, with the value + every candidate carries. A fixed parameter leaves the decision vector, + so the search loses one dimension, and keeps its value in the nine + parameters of every run and of the result. ``None``, the default, + searches every parameter. + :type fixed: dict[str, float] | None + + :param stations: Ids of the stations whose efficiency the objective averages. + ``None``, the default, uses every station the observed series and the + run have in common. + :type stations: tuple[str, ...] | None + + :param allow_blocking_problems: Whether the calibration loads, and searches + on, a configuration whose inputs carry blocking problems. The checks + still run and every problem is still reported; the search then starts + instead of stopping. ``False``, the default, refuses such a + configuration, as a run of the model does. The value is recorded in + :file:`result.json`, since what the search fitted was measured on inputs + the validation rejected. + :type allow_blocking_problems: bool + """ + + variable: str = "arn" + spinup_steps: int = 0 + maxiter: int = 100 + popsize: int = 15 + seed: int | None = None + workers: int | None = None + temp_dir: str | None = None + init: str | np.ndarray = "sobol" + strategy: str = "best1exp" + mutation: tuple[float, float] = (0.5, 1.0) + recombination: float = 0.7 + polish: bool = False + bounds: dict[str, tuple[float, float]] | None = None + fixed: dict[str, float] | None = None + stations: tuple[str, ...] | None = None + allow_blocking_problems: bool = False + + +@dataclass(frozen=True) +class CalibrationResult: + """What a finished calibration found and where it wrote it. + + :param best_parameters: The nine calibration parameters of the best + candidate, the derived ``w_3`` included. + :type best_parameters: dict[str, float] + + :param best_nse: The mean Nash-Sutcliffe efficiency of the best candidate, + ``None`` when no recorded evaluation reached one. + :type best_nse: float | None + + :param best_objective: The objective value of the best candidate. + :type best_objective: float + + :param evaluations: Number of objective evaluations the search spent. + :type evaluations: int + + :param generations: Number of generations the search ran. + :type generations: int + + :param success: Whether SciPy reports the search as successful. + :type success: bool + + :param message: The termination message of SciPy. + :type message: str + + :param run_dir: The run directory the artifacts were written to. + :type run_dir: pathlib.Path + + :param evaluations_csv: The table of every evaluation. + :type evaluations_csv: pathlib.Path + + :param result_json: The summary of the calibration. + :type result_json: pathlib.Path + + :param calibrated_config: The input configuration with the best parameters. + :type calibrated_config: pathlib.Path + + :param best_series: The observed and the simulated series of the best + candidate, side by side on the compared steps. + :type best_series: pathlib.Path + + :param stations_csv: The goodness of fit of the best candidate at every + station the two series share. + :type stations_csv: pathlib.Path + """ + + best_parameters: dict[str, float] + best_nse: float | None + best_objective: float + evaluations: int + generations: int + success: bool + message: str + run_dir: Path + evaluations_csv: Path + result_json: Path + calibrated_config: Path + best_series: Path + stations_csv: Path + + +def _prepare_worker() -> None: + """Silence a worker process before it evaluates its candidate. + + The pool starts every worker as a fresh interpreter with no logging + configuration, so the warnings the loader emits for a configuration whose + inputs are not validated again would reach the terminal through the last + resort handler of :mod:`logging`, once per evaluation; the parent already + reported them when it validated the inputs. The PCRaster framework also + prints a newline when an interpreter exits, one per evaluation. A null + handler on the root logger and a standard output pointed at the null device + keep both out of the terminal: every failure of an evaluation is recorded in + its JSON record and summarized by the parent instead. + """ + logging.getLogger().addHandler(logging.NullHandler()) + sys.stdout = Path(os.devnull).open("w", encoding="utf-8") + + +def _report_failures(records: list[dict[str, Any]]) -> None: + """Log one summary of the evaluations that failed, if any did.""" + failures = [ + record["error"] + for record in records + if record.get("error") not in (None, _INADMISSIBLE_ERROR) + ] + if failures: + logger.warning( + "%d of %d evaluation(s) failed and were ranked behind every evaluated candidate; " + "the first failure was: %s. Every failure is recorded in the evaluations table.", + len(failures), + len(records), + failures[0], + ) + + +@dataclass +class _Progress: + """The per-generation callback: it only logs, it never stops the search.""" + + evaluations_dir: Path + generations: int = field(default=0) + + def __call__(self, intermediate_result) -> bool: + """Log the generation, the best objective and efficiency so far and the records. + + The best efficiency is read from the records the workers have written, + which is where the efficiencies live: what SciPy reports is the + objective, and the efficiency behind it is what a person watching a + calibration reads. + + SciPy inspects the signature of the callback and passes the + intermediate result by keyword under exactly this name; returning + ``True`` would halt the search, so the callback always returns + ``False``. + """ + self.generations = int(getattr(intermediate_result, "nit", 0) or 0) + records = _read_records(self.evaluations_dir) + best = _best_record(records) + best_nse = None if best is None else best.get("nse") + progress.info( + "Generation %d: best objective %.6g, best NSE %s, %d evaluation(s) recorded.", + self.generations, + float(intermediate_result.fun), + "n/a" if best_nse is None else f"{float(best_nse):.6f}", + len(records), + ) + return False + + +def calibrate( + config_path: PathInput, + observed_path: PathInput, + run_dir: PathInput, + settings: CalibrationSettings = CalibrationSettings(), +) -> CalibrationResult: + """Calibrate a configuration against an observed series and write the artifacts. + + :param config_path: The configuration file to calibrate, legacy or format 1.0. + :type config_path: str | os.PathLike[str] + + :param observed_path: The observed series, in the CSV layout of the time + series tables of the model or as a PCRaster TSS file. + :type observed_path: str | os.PathLike[str] + + :param run_dir: Directory the artifacts of the calibration are written to. + It is created when it does not exist. + :type run_dir: str | os.PathLike[str] + + :param settings: How the search is run, defaults to :class:`CalibrationSettings`. + :type settings: CalibrationSettings, optional + + :return: The best candidate and the artifacts of the run. + :rtype: CalibrationResult + + :raises CalibrationError: If SciPy is not installed, if the configuration + cannot be calibrated as asked, if the observed series shares no time + step or no station with what the configuration will sample, if the + station selection names none of those stations, if a worker process + died before finishing its evaluation, or if every evaluation of the run + failed. An interrupted or crashed search still leaves the table of the + evaluations it made behind. + :raises ValueError: If the decision space the settings describe is not a + search (an unknown or derived parameter, a fixed value or an overridden + bound outside the range of the application settings, every parameter + fixed), or if an explicit initial population does not have one column per + free parameter. + :raises ImportError: If PCRaster or GDAL are not installed. + :raises FileNotFoundError: If the configuration or the observed series is + not there. + :raises ConfigurationError: If the configuration carries blocking problems. + """ + differential_evolution = _require_scipy() + + from ..api import Model + + if settings.variable not in VARIABLE_IDS: + raise CalibrationError( + f"'{settings.variable}' is not an output variable of the model; " + f"calibrate against one of {', '.join(VARIABLE_IDS)}." + ) + if settings.spinup_steps < 0: + raise CalibrationError( + f"The number of spin-up steps must not be negative, got {settings.spinup_steps}." + ) + if settings.workers is not None and settings.workers < 1: + raise CalibrationError( + f"The calibration needs at least one worker process, got {settings.workers}." + ) + _check_stations(settings.stations) + + space = decision_space(fixed=settings.fixed, bounds=settings.bounds) + if not isinstance(settings.init, str): + shape = np.shape(settings.init) + if len(shape) != 2 or shape[1] != space.dimension: + raise ValueError( + f"The initial population has the shape {shape} and the search has " + f"{space.dimension} free parameter(s) ({', '.join(space.free_names)}); it must " + f"be a two-dimensional array of {shape[0] if shape else 0} member(s) with one " + "column per free parameter." + ) + + config_file = as_path(config_path).absolute() + # The only validation of a whole calibration: the workers never revalidate, + # so this is where a blocking problem stops the search, and where + # ``allow_blocking_problems`` lets it start anyway. + configuration = Model.from_file( + config_file, + validate_input=True, + allow_blocking_problems=settings.allow_blocking_problems, + ).configuration + file_v1 = _as_v1(configuration) + _check_sample_locations(configuration, file_v1) + + observed = read_series(observed_path) + aggregation = file_v1.model_simulation_output.time_series_samples.aggregation + sampled_ids = ( + None + if aggregation is Aggregation.ZONES + else _station_ids(configuration.raster_files.sample_locations) + ) + compared_steps = _compared_steps( + first_step=configuration.simulation_period.first_step, + last_step=configuration.simulation_period.last_step, + spinup_steps=settings.spinup_steps, + ) + _check_observed_series( + observed, + first_step=configuration.simulation_period.first_step, + last_step=configuration.simulation_period.last_step, + spinup_steps=settings.spinup_steps, + aggregation=aggregation, + station_ids=sampled_ids, + ) + _check_station_selection(settings.stations, observed, sampled_ids) + x0 = _starting_point(space, configuration.calibration_parameters.model_dump()) + + run_directory = as_path(run_dir).absolute() + evaluations_dir = run_directory / EVALUATIONS_DIRNAME + evaluations_dir.mkdir(parents=True, exist_ok=True) + previous = _record_paths(evaluations_dir) + if previous: + raise CalibrationError( + f"The run directory {run_directory} already holds {len(previous)} evaluation " + "record(s) of an earlier calibration; consolidating both runs into one table " + "would mix them. Choose an empty run directory." + ) + observed_csv = run_directory / OBSERVED_CSV + _write_observed_csv(observed_csv, observed, compared_steps, settings.stations) + _report_observed_coverage(observed, compared_steps, settings.stations) + temp_dir = ( + as_path(settings.temp_dir).absolute() + if settings.temp_dir is not None + else Path(tempfile.gettempdir()) + ) + temp_dir.mkdir(parents=True, exist_ok=True) + + workers = ( + settings.workers if settings.workers is not None else max(1, (os.cpu_count() or 2) - 1) + ) + context = EvaluationContext( + document=file_v1.to_dict(), + base_dir=configuration.base_dir, + variable=settings.variable, + observed=observed, + spinup_steps=settings.spinup_steps, + temp_dir=str(temp_dir), + evaluations_dir=str(evaluations_dir), + validate_input=False, + space=space, + stations=settings.stations, + ) + + planned_members = _population_size(settings, space.dimension) + progress.info( + "Calibrating %d free parameter(s) with %d population member(s) per generation and " + "at most %d generation(s): up to %d model run(s), on %d worker process(es).", + space.dimension, + planned_members, + settings.maxiter, + planned_members * (settings.maxiter + 1), + workers, + ) + _report_decision_space(space, settings.stations) + + # A search with at most one free weight has no constraint to hand over: the + # sum of the weights is then a bound of the decision space, and SciPy takes + # a different code path for a constrained search, so the argument is left + # out instead of being passed as ``None``. + constraint = space.weights_constraint() + callback = _Progress(evaluations_dir) + executor = _pool(workers) + arguments: dict[str, Any] = { + "strategy": settings.strategy, + "maxiter": settings.maxiter, + "popsize": settings.popsize, + "mutation": settings.mutation, + "recombination": settings.recombination, + "rng": settings.seed, + "polish": settings.polish, + "init": settings.init, + "x0": x0, + "updating": "deferred", + "workers": executor.map, + "callback": callback, + } + if constraint is not None: + arguments["constraints"] = constraint + evaluations_csv = run_directory / EVALUATIONS_CSV + try: + try: + optimum = differential_evolution( + functools.partial(evaluate, context=context), + list(space.bounds), + **arguments, + ) + finally: + executor.shutdown(wait=True) + except BaseException as error: + # An interruption or a crash after the first evaluation still leaves a + # table: the evaluations that were made are hours of model runs, and a + # calibration that has to be stopped is read from them like any other. + _consolidate(evaluations_csv, evaluations_dir) + if isinstance(error, BrokenProcessPool): + raise CalibrationError( + "A worker process died before finishing its evaluation, most often " + "because the system killed it for lack of memory: every worker holds " + "the rasters of one model run. Run the calibration again with a lower " + f"--workers. The evaluations made so far are in {evaluations_csv}." + ) from error + raise + + records = _read_records(evaluations_dir) + _write_evaluations_csv(evaluations_csv, records) + + _report_failures(records) + best = _best_record(records) + if best is None: + raise CalibrationError( + "Every evaluation of the calibration failed; the first failure was: " + f"{_first_error(records)}. The table of the evaluations is {evaluations_csv}." + ) + + # The size of the population SciPy actually built, not the size this module + # predicted for the budget: a named initialization resizes the population + # (``sobol`` rounds it up to a power of two) and the summary must describe + # the search that ran. + population_size = int(np.shape(optimum.population)[0]) + if population_size != planned_members: + logger.info( + "The search used %d population member(s), %d were announced.", + population_size, + planned_members, + ) + + best_parameters = space.to_parameters(optimum.x) + best_objective = float(optimum.fun) + result_json = run_directory / RESULT_JSON + calibrated_config = run_directory / f"{config_file.stem}-calibrated.json" + stations_csv = run_directory / STATIONS_CSV + best_series = run_directory / f"best_{settings.variable}.csv" + _write_stations_csv(stations_csv, best.get("station_metrics") or {}, settings.stations) + _write_calibrated_config(calibrated_config, configuration, file_v1, best_parameters) + _write_best_series( + best_series, + document=file_v1.to_dict(), + base_dir=configuration.base_dir, + parameters=best_parameters, + variable=settings.variable, + observed=observed, + spinup_steps=settings.spinup_steps, + temp_dir=temp_dir, + workers=workers, + evaluations_csv=evaluations_csv, + ) + _write_result_json( + result_json, + best_parameters=best_parameters, + best_nse=best["nse"], + best_objective=best_objective, + optimum=optimum, + settings=settings, + space=space, + workers=workers, + temp_dir=temp_dir, + population_size=population_size, + artifacts={ + "evaluations_csv": str(evaluations_csv), + "observed_csv": str(observed_csv), + "stations_csv": str(stations_csv), + "best_series": str(best_series), + "calibrated_config": str(calibrated_config), + }, + ) + + progress.info( + "Calibration finished after %d evaluation(s) in %d generation(s): best objective " + "%.6g, best NSE %s.", + int(optimum.nfev), + int(optimum.nit), + best_objective, + "n/a" if best["nse"] is None else f"{float(best['nse']):.6f}", + ) + return CalibrationResult( + best_parameters=best_parameters, + best_nse=best["nse"], + best_objective=best_objective, + evaluations=int(optimum.nfev), + generations=int(optimum.nit), + success=bool(optimum.success), + message=str(optimum.message), + run_dir=run_directory, + evaluations_csv=evaluations_csv, + result_json=result_json, + calibrated_config=calibrated_config, + best_series=best_series, + stations_csv=stations_csv, + ) + + +def _require_scipy(): + """Return ``scipy.optimize.differential_evolution``, or explain how to install it. + + :raises CalibrationError: If SciPy is not installed. + """ + from .. import _deps + + missing = _deps.missing_calibration_deps() + if missing: + raise CalibrationError(_deps.calibration_deps_message(missing)) + from scipy.optimize import differential_evolution + + return differential_evolution + + +def _check_stations(stations: tuple[str, ...] | None) -> None: + """Refuse a station selection the objective could not be averaged over. + + :param stations: The ids the caller asked the objective to average, or + ``None`` for every station the two series share. + :type stations: tuple[str, ...] | None + + :raises CalibrationError: If the selection is not a non-empty sequence of + non-empty station ids. + """ + if stations is None: + return + if isinstance(stations, str) or not isinstance(stations, (tuple, list)): + raise CalibrationError( + "The stations of the objective must be given as a tuple of station ids, " + f"got {type(stations).__name__}." + ) + if not stations: + raise CalibrationError( + "The stations of the objective are an empty selection; name at least one station, " + "or leave them unset to average every station the two series share." + ) + wrong = [station for station in stations if not isinstance(station, str) or not station.strip()] + if wrong: + raise CalibrationError( + f"The station id(s) {', '.join(repr(station) for station in wrong)} are not " + "non-empty station ids; a station is named by the id it carries in the header of " + "the time series table." + ) + + +def _pool(workers: int) -> ProcessPoolExecutor: + """Return the process pool one evaluation at a time is run in. + + ``spawn`` with one task per worker is what isolates the runs from each + other: PCRaster keeps its clone and its raster memory process-wide, so a + worker reused for a second evaluation would carry the state of the first. + """ + return ProcessPoolExecutor( + max_workers=workers, + mp_context=multiprocessing.get_context("spawn"), + max_tasks_per_child=1, + initializer=_prepare_worker, + ) + + +def _consolidate(evaluations_csv: Path, evaluations_dir: Path) -> None: + """Write the table of the evaluations recorded so far, if there are any. + + This is what an interrupted or crashed search leaves behind: the records + are already on disk, one per evaluation, and consolidating them costs + nothing next to the model runs they stand for. A failure to write the table + is reported and swallowed, since the exception that brought the search down + is the one the caller has to see. + """ + try: + records = _read_records(evaluations_dir) + if records: + _write_evaluations_csv(evaluations_csv, records) + logger.warning( + "The search ended early; the %d evaluation(s) it recorded are in %s.", + len(records), + evaluations_csv, + ) + except Exception as error: # pragma: no cover - the original failure must win + logger.warning("The evaluations of the interrupted search could not be written: %s", error) + + +def _compared_steps(*, first_step: int, last_step: int, spinup_steps: int) -> list[int]: + """Return the simulated time steps the efficiency is computed on, in order.""" + return [step for step in range(first_step, last_step + 1) if step > spinup_steps] + + +def _check_station_selection( + stations: tuple[str, ...] | None, + observed: Series, + sampled_ids: set[str] | None, +) -> None: + """Refuse a station selection that names none of the stations of the run. + + The selection is applied inside every evaluation, so a misspelt id would + make every one of them fail and the whole budget would be spent finding + out. The comparable stations are the ones the observed series and the run + have in common; under the ``zones`` aggregation the parent does not know + the ids of the run, so the selection is only checked against the + observations. + + :param stations: The ids the objective averages, or ``None``. + :type stations: tuple[str, ...] | None + + :param observed: The observed series, as the parent read it. + :type observed: rubem.calibration.objective.Series + + :param sampled_ids: The ids the run will sample, or ``None`` when they are + not known before the run. + :type sampled_ids: set[str] | None + + :raises CalibrationError: If none of the selected ids is comparable. + """ + if stations is None: + return + comparable = set(observed.stations) + if sampled_ids is not None: + comparable &= sampled_ids + unknown = [station for station in stations if station not in comparable] + if len(unknown) == len(stations): + raise CalibrationError( + f"The station(s) {', '.join(stations)} of the objective are none of the " + f"station(s) {', '.join(sorted(comparable)) or 'none'} the observed series and " + "the run have in common; every evaluation would fail." + ) + if unknown: + logger.warning( + "The objective names %d station(s) the comparison does not have (%s); it is " + "averaged over the remaining ones.", + len(unknown), + ", ".join(unknown), + ) + + +def _starting_point(space: DecisionSpace, parameters: dict[str, float]) -> np.ndarray: + """Return the ``x0`` of the search: the configuration, inside the bounds of the run. + + The starting point is the parameter set of the configuration being + calibrated. A bound narrowed for this run may exclude it, and SciPy refuses + a starting point outside the bounds, so every coordinate that falls outside + is moved onto the bound it crosses and the move is reported: the run then + starts from the admissible point closest to the configuration instead of + ending before its first evaluation. + + :param space: The decision vector of the run. + :type space: rubem.calibration.parameters.DecisionSpace + + :param parameters: The calibration parameters of the configuration. + :type parameters: dict[str, float] + + :return: The free values of the starting point, in the order of the space. + :rtype: numpy.ndarray + """ + x0 = space.from_parameters(parameters) + for index, (name, (minimum, maximum)) in enumerate( + zip(space.free_names, space.bounds, strict=True) + ): + value = float(x0[index]) + moved = min(max(value, minimum), maximum) + if moved != value: + x0[index] = moved + logger.warning( + "The configuration has %s = %g, outside the bound [%g, %g] of this " + "calibration; the search starts from %g instead.", + name, + value, + minimum, + maximum, + moved, + ) + if not space.is_admissible(x0): + logger.warning( + "The parameters of the configuration are not an admissible candidate of this " + "calibration; the search starts from the rest of the initial population." + ) + return x0 + + +def _report_decision_space(space: DecisionSpace, stations: tuple[str, ...] | None) -> None: + """Log what this run searches: the fixed parameters, the narrowed bounds and the stations.""" + if space.fixed: + logger.info( + "Fixed, not searched: %s.", + ", ".join(f"{name} = {value:g}" for name, value in space.fixed.items()), + ) + narrowed = [ + (name, bound) + for name, bound in zip(space.free_names, space.bounds, strict=True) + if bound != variable_range(name) + ] + if narrowed: + logger.info( + "Searched in a narrowed range: %s.", + ", ".join( + f"{name} in [{minimum:g}, {maximum:g}]" for name, (minimum, maximum) in narrowed + ), + ) + if stations is not None: + logger.info("The objective averages the station(s) %s.", ", ".join(stations)) + + +def _as_v1(configuration) -> ModelConfigurationFileV1: + """Return the configuration as a format 1.0 model, converting a legacy one.""" + if configuration.file_v1 is not None: + return configuration.file_v1 + return ModelConfigurationFileV1.from_legacy(configuration.file) + + +def _check_sample_locations(configuration, file_v1: ModelConfigurationFileV1) -> None: + """Refuse a configuration whose time series would have nowhere to be sampled. + + :raises CalibrationError: If the raster the configured aggregation needs is + not part of the configuration. + """ + aggregation = file_v1.model_simulation_output.time_series_samples.aggregation + if aggregation is Aggregation.ZONES: + if not configuration.raster_files.zones: + raise CalibrationError( + "The configuration aggregates its time series over zones but names no " + "zones raster; the calibration has no station to compare." + ) + return + if not configuration.raster_files.sample_locations: + raise CalibrationError( + "The configuration names no sample locations raster; the calibration has no " + "station to compare the observed series with." + ) + + +def _station_ids(sample_locations: PathInput) -> set[str]: + """Return the station ids of a sample locations raster. + + The raster is a nominal map, so its cells read as integers. Neither the + no-data cells nor the zeroes are stations, which is what the validation of + the inputs already asks of a sample raster: its ids are ``1..N`` and its + background is missing or zero. Every remaining value becomes the string the + model writes in the header of its time series table, ``str(int(value))``. + + :param sample_locations: The sample locations raster of the configuration. + :type sample_locations: str | os.PathLike[str] + + :return: The distinct station ids, as strings. + :rtype: set[str] + """ + from osgeo import gdal + + gdal.UseExceptions() + gdal.AllRegister() + dataset = gdal.Open(str(as_path(sample_locations))) + try: + band = dataset.GetRasterBand(1) + values = np.unique(np.asarray(band.ReadAsArray())) + no_data = band.GetNoDataValue() + finally: + dataset = None + + if no_data is not None: + values = values[values != no_data] + return {str(int(value)) for value in values.tolist() if math.isfinite(value) and value != 0} + + +def _check_observed_series( + observed: Series, + *, + first_step: int, + last_step: int, + spinup_steps: int, + aggregation: Aggregation, + station_ids: set[str] | None, +) -> None: + """Refuse an observed series that cannot be compared with what the run samples. + + Every evaluation aligns the observed series with the table its own + simulation wrote, so a series that shares no step or no station with that + table makes every evaluation fail. The mismatch is the same for every + candidate, so it is settled here, in the parent, before a single model run + is spent on it. + + The compared steps are the simulated steps, ``first_step`` to ``last_step``, + that survive the spin-up. The stations are the ids the ``point`` and the + ``subcatchment`` aggregations sample, the values of the sample locations + raster; the ``zones`` aggregation remaps the zone ids to ``1..N`` at run + time (the mapping is written to ``zones_mapping.csv``), so the parent cannot + know the ids of its stations and the station check is skipped for it. + + :param observed: The observed series, as the parent read it. + :type observed: rubem.calibration.objective.Series + + :param first_step: Number of the first simulated step. + :type first_step: int + + :param last_step: Number of the last simulated step. + :type last_step: int + + :param spinup_steps: Number of initial time steps excluded from the + efficiency. + :type spinup_steps: int + + :param aggregation: The aggregation the time series of the run use. + :type aggregation: Aggregation + + :param station_ids: The ids of the stations the run will sample, or ``None`` + when they are not known before the run. + :type station_ids: set[str] | None + + :raises CalibrationError: If no observed step is compared, or if no observed + station id is one of the stations of the run. + """ + compared_steps = {step for step in range(first_step, last_step + 1) if step > spinup_steps} + observed_steps = {int(step) for step in np.asarray(observed.steps).tolist()} + if not compared_steps & observed_steps: + raise CalibrationError( + f"The observed series covers the step(s) {min(observed_steps)} to " + f"{max(observed_steps)}, and the configuration simulates the step(s) " + f"{first_step} to {last_step}, of which the {spinup_steps} spin-up step(s) " + "are excluded from the efficiency: the two series would have no time step " + "in common and every evaluation would fail." + ) + + if aggregation is Aggregation.ZONES or station_ids is None: + return + + observed_ids = set(observed.stations) + if not observed_ids & station_ids: + raise CalibrationError( + f"The observed series has the station(s) {', '.join(sorted(observed_ids)) or 'none'} " + f"and the configuration samples the station(s) " + f"{', '.join(sorted(station_ids)) or 'none'}: the two series would have no " + "station in common and every evaluation would fail." + ) + + foreign = sorted(observed_ids - station_ids) + if foreign: + logger.warning( + "The observed series has %d station(s) the configuration does not sample " + "(%s); they are ignored by the calibration.", + len(foreign), + ", ".join(foreign), + ) + + +def _observed_window(observed: Series, compared_steps: list[int]) -> np.ndarray: + """Return the rows of the observed series that fall on the compared steps.""" + steps = np.asarray(observed.steps, dtype=np.int64) + return np.flatnonzero(np.isin(steps, np.asarray(compared_steps, dtype=np.int64))) + + +def _in_selection(station: str, stations: tuple[str, ...] | None) -> str: + """Return whether a station enters the objective, as the tables spell it.""" + return "true" if stations is None or station in stations else "false" + + +def _write_observed_csv( + path: Path, + observed: Series, + compared_steps: list[int], + stations: tuple[str, ...] | None, +) -> None: + """Write what the observations offer on the compared steps, one row per station. + + The row of a station carries how many of the compared steps it observes + (``pairs_in_window``), how many values on those steps the mask of + :func:`rubem.calibration.objective.valid_mask` rejected as gaps + (``dropped``, the negative values among them) and the summary statistics of + the values that were kept. It is written before the search: a station whose + gauge is nearly empty over the simulated window explains a mean efficiency + afterwards, and it is cheaper to read it first. + """ + window = _observed_window(observed, compared_steps) + with path.open("w", encoding="utf-8", newline="") as table: + writer = csv.writer(table, delimiter=STATION_SEPARATOR) + writer.writerow(OBSERVED_COLUMNS) + for station in sorted(observed.stations): + values = np.asarray(observed.stations[station], dtype=np.float64)[window] + mask = valid_mask(values) + kept = values[mask] + writer.writerow( + [ + station, + _in_selection(station, stations), + int(kept.size), + int(values.size - kept.size), + _cell(float(np.mean(kept)) if kept.size else None), + _cell(float(np.std(kept, ddof=1)) if kept.size > 1 else None), + _cell(float(np.min(kept)) if kept.size else None), + _cell(float(np.max(kept)) if kept.size else None), + ] + ) + + +def _report_observed_coverage( + observed: Series, + compared_steps: list[int], + stations: tuple[str, ...] | None, +) -> None: + """Log how much of the compared window the observations cover, per station.""" + window = _observed_window(observed, compared_steps) + covered = np.asarray(observed.steps, dtype=np.int64)[window] + if covered.size: + span = f"the step(s) {int(covered.min())} to {int(covered.max())}" + else: + span = "no step of it" + progress.info( + "The calibration compares %d time step(s); the observed series covers %d of them, %s.", + len(compared_steps), + int(covered.size), + span, + ) + for station in sorted(observed.stations): + values = np.asarray(observed.stations[station], dtype=np.float64)[window] + kept = int(np.count_nonzero(valid_mask(values))) + progress.info( + "Station %s%s: %d observed value(s) on the compared steps, %d dropped as gaps.", + station, + "" if stations is None or station in stations else " (not in the objective)", + kept, + int(values.size) - kept, + ) + + +def _write_stations_csv( + path: Path, + station_metrics: dict[str, Any], + stations: tuple[str, ...] | None, +) -> None: + """Write the goodness of fit of the best candidate, one row per station. + + The statistics are the ones its own evaluation recorded, not a + recomputation: the table and the record of that evaluation therefore always + agree. A station outside the selection is measured like any other and + marked as such, which is what makes the stations left out of the objective a + validation of the calibration. + """ + with path.open("w", encoding="utf-8", newline="") as table: + writer = csv.writer(table, delimiter=STATION_SEPARATOR) + writer.writerow(STATION_COLUMNS) + for station in sorted(station_metrics): + measured = station_metrics[station] or {} + writer.writerow( + [ + station, + _in_selection(station, stations), + int(measured.get("pairs") or 0), + # The remaining columns are the statistics, under the very + # names the record spells them with. + *(_cell(measured.get(name)) for name in STATION_COLUMNS[3:]), + ] + ) + + +def _write_best_series( + path: Path, + *, + document: dict, + base_dir: str | None, + parameters: dict[str, float], + variable: str, + observed: Series, + spinup_steps: int, + temp_dir: Path, + workers: int, + evaluations_csv: Path, +) -> None: + """Run the best candidate once more and write its series beside the observed one. + + The search keeps no output: every evaluation writes into a temporary + directory that is removed when it ends, which is what keeps a calibration of + thousands of runs from filling a disk. The candidate that won is therefore + run one last time, in a worker of its own, so that the series behind the + reported efficiency can be plotted against the observations. The steps and + the stations of the table are the ones the objective compared, from the same + alignment, and a value either series does not offer is left as an empty + cell. + + :raises CalibrationError: If the run of the best candidate failed. + """ + output_dir = Path(tempfile.mkdtemp(prefix="best-", dir=str(temp_dir))) + executor = _pool(workers) + try: + table = executor.submit( + simulate_best, document, base_dir, parameters, str(output_dir), variable + ).result() + simulated = read_series(table) + except Exception as error: + raise CalibrationError( + "The best candidate of the search could not be run again to keep its series: " + f"{type(error).__name__}: {error}. The calibration itself finished; its " + f"evaluations are in {evaluations_csv}." + ) from error + finally: + executor.shutdown(wait=True) + shutil.rmtree(output_dir, ignore_errors=True) + + alignment = align_series(simulated, observed, spinup_steps) + columns = [ + column + for station in alignment.stations + for column in (f"observed_{station}", f"simulated_{station}") + ] + with path.open("w", encoding="utf-8", newline="") as best: + writer = csv.writer(best, delimiter=STATION_SEPARATOR) + writer.writerow(["step", *columns]) + for row, step in enumerate(alignment.steps.tolist()): + cells: list[str] = [] + for station in alignment.stations: + cells.append( + _observation(observed.stations[station], alignment.observed_index[row]) + ) + cells.append( + _observation(simulated.stations[station], alignment.simulated_index[row]) + ) + writer.writerow([int(step), *cells]) + + +def _observation(values: np.ndarray, row: int) -> str: + """Return one value of a series, or an empty cell when the mask rejects it.""" + value = float(np.asarray(values, dtype=np.float64)[row]) + return _cell(value) if bool(valid_mask(np.asarray([value]))[0]) else "" + + +def _population_size(settings: CalibrationSettings, dimension: int = len(FREE_PARAMETERS)) -> int: + """Return the number of population members the search will use. + + The arithmetic mirrors SciPy's: an explicit population is used as it is, a + named initialization gives ``max(5, popsize * dimension)`` members, with + ``dimension`` the number of free parameters, and ``sobol`` alone rounds that + up to the next power of two. + """ + if not isinstance(settings.init, str): + return int(np.shape(settings.init)[0]) + members = max(5, settings.popsize * dimension) + if settings.init == "sobol": + return 1 << (members - 1).bit_length() + return members + + +def _record_paths(evaluations_dir: Path) -> list[Path]: + """Return the JSON records written so far, without the files still being staged.""" + return sorted(path for path in evaluations_dir.glob("*.json") if not path.name.startswith(".")) + + +def _read_records(evaluations_dir: Path) -> list[dict[str, Any]]: + """Read every JSON record of the run, in the order of :func:`_record_order`. + + A record that cannot be read is reported and skipped: an unreadable record + of one evaluation must not hide the result of all the others. + """ + records = [] + for path in _record_paths(evaluations_dir): + try: + records.append(json.loads(path.read_text(encoding="utf-8"))) + except (OSError, ValueError) as error: + logger.warning("Skipping the unreadable evaluation record %s: %s", path, error) + records.sort(key=_record_order) + return records + + +def _record_order(record: dict[str, Any]) -> tuple[tuple[float, ...], str]: + """Return the sort key of a record: its parameters, then its identifier. + + The records are not ordered by the moment they were written: the workers run + in parallel, so that order changes from run to run. Ordering them by the + candidate they evaluated takes that scheduling order out of + ``evaluations.csv``. It does not make the table reproducible on its own: the + ``id``, ``pid`` and ``elapsed_seconds`` columns differ between two runs, and + from the first generation on the selection reads the objective values, so + two runs of one seed evaluate the same candidates only when the simulation + itself is reproducible (on the real basins that requires a fixed + ``RASTERS.ldd``, since ``lddcreate`` does not always derive the same + directions twice; see the note on the LDD raster in the user guide). The + identifier only separates two evaluations of the very same candidate. + """ + parameters = record.get("parameters") or {} + values = tuple(float(parameters.get(name, math.inf)) for name in CALIBRATION_PARAMETERS) + return values, str(record.get("id", "")) + + +def _write_evaluations_csv(path: Path, records: list[dict[str, Any]]) -> None: + """Write one row per evaluation, in the columns of :data:`EVALUATION_COLUMNS`.""" + with path.open("w", encoding="utf-8", newline="") as table: + writer = csv.writer(table) + writer.writerow(EVALUATION_COLUMNS) + for record in records: + parameters = record.get("parameters") or {} + writer.writerow( + [ + record.get("id", ""), + record.get("pid", ""), + record.get("started_at", ""), + *(_cell(parameters.get(name)) for name in CALIBRATION_PARAMETERS), + _cell(record.get("nse")), + _cell(record.get("objective")), + _cell(record.get("elapsed_seconds")), + record.get("error") or "", + ] + ) + + +def _cell(value) -> str: + """Return a numeric cell with full precision, or an empty cell for a missing value.""" + return "" if value is None else repr(float(value)) + + +def _best_record(records: list[dict[str, Any]]) -> dict[str, Any] | None: + """Return the record of the lowest objective among those that were evaluated. + + A record without an objective is not a candidate for the best one: the + records are written by other processes, and the parent reads whatever is in + the directory rather than assuming every file there is complete. + """ + evaluated = [ + record + for record in records + if not record.get("error") and record.get("objective") is not None + ] + if not evaluated: + return None + return min(evaluated, key=lambda record: float(record["objective"])) + + +def _first_error(records: list[dict[str, Any]]) -> str: + """Return the error of the first record that carries one.""" + return next((str(record["error"]) for record in records if record.get("error")), "unknown") + + +def _write_result_json( + path: Path, + *, + best_parameters: dict[str, float], + best_nse: float | None, + best_objective: float, + optimum, + settings: CalibrationSettings, + space: DecisionSpace, + workers: int, + temp_dir: Path, + population_size: int, + artifacts: dict[str, str], +) -> None: + """Write the summary of the calibration. + + The ``artifacts`` object names every file the run wrote next to this one, so + that the summary is the one file a reader has to open to find the rest. + """ + document = { + "best_parameters": best_parameters, + "best_nse": best_nse, + "best_objective": best_objective, + "nfev": int(optimum.nfev), + "nit": int(optimum.nit), + "success": bool(optimum.success), + "message": str(optimum.message), + "population_size": population_size, + "settings": _settings_document(settings, space, workers, temp_dir), + "artifacts": artifacts, + } + path.write_text(json.dumps(document, indent=2) + "\n", encoding="utf-8") + + +def _settings_document( + settings: CalibrationSettings, space: DecisionSpace, workers: int, temp_dir: Path +) -> dict[str, Any]: + """Return the settings as JSON-able data, resolved the way the run used them. + + The number of workers, the temporary directory and the bounds are the values + the run actually used, not the ``None`` the caller may have left them at, so + that the summary describes a calibration that can be repeated: the bounds are + the effective ones, one entry per searched parameter, the ranges of the + application settings narrowed by whatever the caller asked for. + """ + init = settings.init if isinstance(settings.init, str) else np.asarray(settings.init).tolist() + return { + "variable": settings.variable, + "spinup_steps": settings.spinup_steps, + "maxiter": settings.maxiter, + "popsize": settings.popsize, + "seed": settings.seed, + "workers": workers, + "temp_dir": str(temp_dir), + "init": init, + "strategy": settings.strategy, + "mutation": list(settings.mutation), + "recombination": settings.recombination, + "polish": settings.polish, + "bounds": { + name: [minimum, maximum] + for name, (minimum, maximum) in zip(space.free_names, space.bounds, strict=True) + }, + "fixed": dict(space.fixed), + "stations": list(settings.stations) if settings.stations is not None else None, + "allow_blocking_problems": settings.allow_blocking_problems, + } + + +def _write_calibrated_config( + path: Path, + configuration, + file_v1: ModelConfigurationFileV1, + best_parameters: dict[str, float], +) -> None: + """Write the calibrated configuration in the format the input was written in. + + The document goes through the configuration models, so what is written is a + configuration the loader accepts. Its paths are the ones the loader resolved, + which are absolute even when the input file wrote them relative to its own + directory. + """ + calibrated = file_v1.model_copy( + update={ + "model_calibration_parameters": type(file_v1.model_calibration_parameters)( + alpha=best_parameters["alpha"], + b=best_parameters["beta"], + w_1=best_parameters["w_1"], + w_2=best_parameters["w_2"], + w_3=best_parameters["w_3"], + rcd=best_parameters["rcd"], + f=best_parameters["f"], + alpha_gw=best_parameters["alpha_gw"], + x=best_parameters["x"], + ) + } + ) + document = ( + calibrated.to_dict() + if configuration.file_v1 is not None + else calibrated.to_legacy().to_dict() + ) + path.write_text(json.dumps(document, indent=2) + "\n", encoding="utf-8") diff --git a/rubem/cli.py b/rubem/cli.py index dcf83d15..cf19f25f 100644 --- a/rubem/cli.py +++ b/rubem/cli.py @@ -1,7 +1,8 @@ """The ``rubem`` command line. ``rubem run -c [-s | --allow-blocking-problems]`` runs a -simulation; ``rubem config schema`` prints the JSON Schema of the configuration +simulation; ``rubem calibrate`` fits the calibration parameters to an observed +series; ``rubem config schema`` prints the JSON Schema of the configuration file. The former ``rubem -c `` spelling still works for one minor release and emits a ``DeprecationWarning``. """ @@ -10,6 +11,7 @@ import logging import logging.config import logging.handlers +import math import sys import warnings from collections.abc import Sequence @@ -21,7 +23,7 @@ from pydantic import ValidationError from . import __release__ -from ._deps import require_runtime_deps +from ._deps import require_calibration_deps, require_runtime_deps from .configuration._problems import ConfigurationError from .configuration.app_settings import AppSettings from .preprocessing.cli import app as preprocess_app @@ -70,6 +72,122 @@ class SchemaFormat(StrEnum): legacy = "legacy" +class InitMethod(StrEnum): + """Initial populations the differential evolution can be started from.""" + + sobol = "sobol" + latinhypercube = "latinhypercube" + halton = "halton" + random = "random" + + +def _parse_bounds(values: list[str] | None) -> dict[str, tuple[float, float]]: + """Return the ``--bound NAME=MIN:MAX`` options as a mapping of ranges. + + Only the shape of the text is checked here; whether the range narrows the + range of the application settings is settled by the calibration itself, + which knows the parameters. + + :param values: The options as they were written, or ``None``. + :type values: list[str] | None + + :return: The ``(minimum, maximum)`` range of each named parameter. + :rtype: dict[str, tuple[float, float]] + + :raises typer.BadParameter: If an option is not ``NAME=MIN:MAX``, or if a + parameter is bounded twice. + """ + parsed: dict[str, tuple[float, float]] = {} + for entry in values or []: + name, separator, rest = entry.partition("=") + minimum, colon, maximum = rest.partition(":") + if not separator or not colon or not name.strip(): + raise typer.BadParameter( + f"'{entry}' is not a bound; write it as NAME=MIN:MAX, as in 'rcd=2:5'.", + param_hint="'--bound'", + ) + parsed[_unique(name.strip(), parsed, "--bound")] = ( + _as_number(minimum, entry, "--bound"), + _as_number(maximum, entry, "--bound"), + ) + return parsed + + +def _parse_fixed(values: list[str] | None) -> dict[str, float]: + """Return the ``--fix NAME=VALUE`` options as a mapping of pinned values. + + :param values: The options as they were written, or ``None``. + :type values: list[str] | None + + :return: The value each named parameter is pinned to. + :rtype: dict[str, float] + + :raises typer.BadParameter: If an option is not ``NAME=VALUE``, or if a + parameter is fixed twice. + """ + parsed: dict[str, float] = {} + for entry in values or []: + name, separator, value = entry.partition("=") + if not separator or not name.strip(): + raise typer.BadParameter( + f"'{entry}' is not a fixed parameter; write it as NAME=VALUE, as in 'x=0'.", + param_hint="'--fix'", + ) + parsed[_unique(name.strip(), parsed, "--fix")] = _as_number(value, entry, "--fix") + return parsed + + +def _parse_stations(value: str | None) -> tuple[str, ...]: + """Return the ``--stations ID[,ID...]`` option as the ids it names. + + :param value: The option as it was written, or ``None``. + :type value: str | None + + :return: The station ids, in the order they were given. + :rtype: tuple[str, ...] + + :raises typer.BadParameter: If the option names no station, or leaves one of + the comma-separated ids empty. + """ + if value is None: + return () + stations = tuple(station.strip() for station in value.split(",")) + if not all(stations): + raise typer.BadParameter( + f"'{value}' does not name the stations; write them as ID[,ID...], as in '1,2,3'.", + param_hint="'--stations'", + ) + return stations + + +def _as_number(text: str, entry: str, option: str) -> float: + """Return a number written on the command line. + + :raises typer.BadParameter: If the text is not a finite number. + """ + try: + value = float(text) + except ValueError as error: + raise typer.BadParameter( + f"'{text.strip()}' of '{entry}' is not a number.", param_hint=f"'{option}'" + ) from error + if not math.isfinite(value): + raise typer.BadParameter( + f"'{text.strip()}' of '{entry}' is not a finite number.", param_hint=f"'{option}'" + ) + return value + + +def _unique(name: str, parsed: dict, option: str) -> str: + """Return a parameter name that was not given twice on the command line. + + :raises typer.BadParameter: If the name is already in the mapping. + """ + if name in parsed: + raise typer.BadParameter(f"'{name}' is given more than once.", param_hint=f"'{option}'") + return name + + def _version_callback(value: bool) -> None: if value: print(f"RUBEM v{__release__}") @@ -197,6 +315,167 @@ def run( logger.info("RUBEM successfully finished!") +@app.command() +def calibrate( + configfile: Annotated[ + Path, + typer.Option( + "-c", + "--configfile", + callback=_configfile_callback, + help="Path to the configuration file (JSON).", + ), + ], + observed: Annotated[ + Path, + typer.Option( + "--observed", + exists=True, + dir_okay=False, + help="Observed series at the sample stations (CSV or PCRaster TSS).", + ), + ], + run_dir: Annotated[ + Path, + typer.Option("-o", "--run-dir", help="Directory for the calibration artifacts."), + ], + variable: Annotated[ + str, typer.Option("--variable", help="Output variable compared with the observed series.") + ] = "arn", + spinup_steps: Annotated[ + int, + typer.Option("--spinup-steps", min=0, help="Leading time steps excluded from the NSE."), + ] = 0, + maxiter: Annotated[ + int, typer.Option("--maxiter", min=1, help="Maximum number of generations.") + ] = 100, + popsize: Annotated[ + int, typer.Option("--popsize", min=1, help="Population multiplier of the search.") + ] = 15, + seed: Annotated[ + int | None, typer.Option("--seed", help="Seed of the differential evolution.") + ] = None, + workers: Annotated[ + int | None, + typer.Option("--workers", min=1, help="Parallel evaluations (default: all cores but one)."), + ] = None, + temp_dir: Annotated[ + Path | None, + typer.Option("--temp-dir", help="Parent of the per-evaluation output directories."), + ] = None, + bound: Annotated[ + list[str] | None, + typer.Option("--bound", help="Narrow the range of one parameter: NAME=MIN:MAX."), + ] = None, + fix: Annotated[ + list[str] | None, + typer.Option("--fix", help="Keep one parameter out of the search: NAME=VALUE."), + ] = None, + stations: Annotated[ + str | None, + typer.Option( + "--stations", + help="Station ids of the objective, comma-separated (default: every shared station).", + ), + ] = None, + init: Annotated[ + InitMethod, typer.Option("--init", help="Initial population of the search.") + ] = InitMethod.sobol, + strategy: Annotated[ + str, typer.Option("--strategy", help="Differential evolution strategy.") + ] = "best1exp", + polish: Annotated[ + bool, + typer.Option("--polish/--no-polish", help="Refine the best candidate with a local search."), + ] = False, + allow_blocking_problems: Annotated[ + bool, + typer.Option( + "--allow-blocking-problems", + help=( + "Search even if the input validation finds blocking problems; " + "they are logged as errors instead of stopping the calibration." + ), + ), + ] = False, +) -> None: + """Calibrate the model parameters against an observed series.""" + # The options that carry a small language of their own are read first, so + # that a typo costs nothing but the message naming the form expected. + bounds = _parse_bounds(bound) + fixed = _parse_fixed(fix) + selected_stations = _parse_stations(stations) + + require_runtime_deps() + require_calibration_deps() + + from .calibration.runner import CalibrationError, CalibrationSettings + from .calibration.runner import calibrate as run_calibration + + settings_arguments: dict = { + "variable": variable, + "spinup_steps": spinup_steps, + "maxiter": maxiter, + "popsize": popsize, + "seed": seed, + "init": init.value, + "strategy": strategy, + "polish": polish, + "allow_blocking_problems": allow_blocking_problems, + } + if workers is not None: + settings_arguments["workers"] = workers + if temp_dir is not None: + settings_arguments["temp_dir"] = str(temp_dir) + if bounds: + settings_arguments["bounds"] = bounds + if fixed: + settings_arguments["fixed"] = fixed + if selected_stations: + settings_arguments["stations"] = selected_stations + + try: + # The library only logs; the progress a command-line calibration is + # expected to show is written here, so that an embedded calibration + # stays silent unless its host configures logging. + print("Loading configuration and validating inputs...", flush=True) + print("Calibration started...", flush=True) + result = run_calibration( + configfile, observed, run_dir, CalibrationSettings(**settings_arguments) + ) + except (ConfigurationError, ValidationError, ValueError) as e: + # A configuration the user can fix: no traceback, the message says what. + logger.critical("Invalid configuration: %s", e) + raise typer.Exit(code=1) from e + except CalibrationError as e: + # The calibration itself refused to run or produced nothing usable; + # the message says what, so no traceback either. + logger.critical("Calibration failed: %s", e) + raise typer.Exit(code=1) from e + except KeyboardInterrupt as e: + logger.critical("RUBEM was interrupted by the user.") + raise typer.Exit(code=2) from e + except Exception as e: + logger.critical("RUBEM unexpectedly quit.") + logger.exception(e) + raise typer.Exit(code=1) from e + + best_nse = "n/a" if result.best_nse is None else f"{result.best_nse:.6f}" + print("Calibration finished successfully!", flush=True) + print(f"Best NSE: {best_nse}", flush=True) + print(f"Best objective: {result.best_objective:.6g}", flush=True) + print(f"Evaluations: {result.evaluations}", flush=True) + print(f"Generations: {result.generations}", flush=True) + print(f"Evaluations table: {result.evaluations_csv}", flush=True) + print(f"Result: {result.result_json}", flush=True) + print(f"Calibrated configuration: {result.calibrated_config}", flush=True) + print(f"Observed summary: {result.run_dir / 'observed.csv'}", flush=True) + print(f"Stations table: {result.stations_csv}", flush=True) + print(f"Best candidate series: {result.best_series}", flush=True) + + logger.info("RUBEM successfully finished!") + + @config_app.command("schema") def config_schema( schema_format: Annotated[ @@ -277,7 +556,7 @@ def _legacy_arguments(argv: Sequence[str]) -> list[str]: global _legacy_invocation _legacy_invocation = False arguments = list(argv) - if not arguments or arguments[0] in ("run", "config"): + if not arguments or arguments[0] in ("run", "config", "calibrate", "preprocess"): return arguments options = [a for a in arguments if a.startswith("-")] if not options or not {o.split("=", 1)[0] for o in options} <= _LEGACY_OPTIONS: diff --git a/tests/conftest.py b/tests/conftest.py index eb861b6f..cd305f0f 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,8 +1,68 @@ import logging import os +from pathlib import Path import pytest +DATASET_DIR_VARIABLE = "RUBEM_DATASET_DIR" +"""Environment variable pointing at the local copy of the published datasets. + +The directory it names holds the manifest ``_basin.json`` of every +published basin beside the extracted ``_basin/`` directory of its inputs. +The datasets are hundreds of megabytes and are not part of the repository, so +the tests that need them are marked ``dataset`` and run only where the variable +is set; everywhere else the suite stays self-contained. +""" + +NO_DATASET_DIR = f"{DATASET_DIR_VARIABLE} is not set" +"""The reason a test marked ``dataset`` is skipped with when the variable is unset.""" + + +def _dataset_dir() -> tuple[Path | None, str]: + """Return the datasets directory, or ``None`` and the reason it is not available. + + An unset or empty variable means the published datasets are not on this + machine; a variable that names something other than a directory is reported + as such, so that a mistyped path or a mount that is not up does not read as + "not requested". + """ + value = os.environ.get(DATASET_DIR_VARIABLE) + if not value: + return None, NO_DATASET_DIR + directory = Path(value) + if not directory.is_dir(): + return None, f"{DATASET_DIR_VARIABLE}={value} is not a directory" + return directory, "" + + +def pytest_collection_modifyitems(config, items): + """Skip every test marked ``dataset`` when the datasets are not there. + + The marker alone is enough: a test that never asks for the ``dataset_dir`` + fixture, because it reaches the datasets through another fixture or through + a helper, is skipped for the same reason as one that does. + """ + directory, reason = _dataset_dir() + if directory is not None: + return + skip = pytest.mark.skip(reason=reason) + for item in items: + if item.get_closest_marker("dataset") is not None: + item.add_marker(skip) + + +@pytest.fixture +def dataset_dir(): + """The directory of the published datasets, or a skipped test. + + :return: The directory :data:`DATASET_DIR_VARIABLE` names. + :rtype: pathlib.Path + """ + directory, reason = _dataset_dir() + if directory is None: + pytest.skip(reason) + return directory + @pytest.fixture(autouse=True) def preserve_working_directory(): diff --git a/tests/dataset/__init__.py b/tests/dataset/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/dataset/test_calibrate_ipojuca.py b/tests/dataset/test_calibrate_ipojuca.py new file mode 100644 index 00000000..4f39c89c --- /dev/null +++ b/tests/dataset/test_calibrate_ipojuca.py @@ -0,0 +1,200 @@ +"""Smoke test of the calibration on the published Ipojuca basin dataset. + +The test rebuilds the configuration of the basin from its manifest, runs the +model once over a short window to obtain the series a perfect calibration has +to reproduce, and then asks the differential evolution to find the parameters +that produced it again, from a five-member initial population, in one +generation. + +The drainage network is fixed before anything runs. ``lddcreate`` is called by +the model whenever the configuration names no ``ldd`` raster, and on the real +basin digital elevation models it returns a different network on every call, so +the accumulated runoff, which is routed along it, is not reproducible between +two runs of one configuration. A calibration on ``arn`` without a fixed network +would compare series routed on different networks; the test therefore generates +the network once and names it in the configuration, as the documentation of the +command requires. + +The dataset is not part of the repository: the whole module is marked +``dataset`` and runs only where ``RUBEM_DATASET_DIR`` points at the local copy. +""" + +import csv +import json +import multiprocessing +import shutil +from pathlib import Path + +import numpy as np +import pytest + +from rubem.api import Model +from rubem.calibration.parameters import CALIBRATION_PARAMETERS, is_admissible +from rubem.calibration.runner import EVALUATION_COLUMNS, CalibrationSettings, calibrate +from rubem.configuration.model_configuration import ModelConfiguration + +# SciPy is the optional ``rubem[calibration]`` extra; an environment without it +# collects this file and skips it instead of failing to import. +pytest.importorskip("scipy.optimize") + +pytestmark = [pytest.mark.dataset, pytest.mark.slow] + +BASIN = "ipojuca_basin" +"""The basin this test calibrates: its manifest is ``.json``.""" + +PATH_SEGMENT = f"{BASIN}/" +"""The segment every input path of the manifest carries. + +The manifests are written by the machine that extracted the datasets and hold +absolute paths of that machine, so every path is rebased onto the local copy by +replacing whatever precedes this segment. +""" + +STEPS = 6 +"""Monthly steps of the window the test simulates, from 01/01/2000 to 01/06/2000.""" + +# Five admissible members around the parameter set of the manifest, the second +# one far from it. SciPy replaces the first member with x0, so the configuration +# under calibration is always part of the initial population, and the search has +# one generation to find it again. +INIT = np.array( + [ + [4.415, 0.078, 0.51, 0.12, 5.375, 0.581, 0.922, 0.307], + [9.500, 0.900, 0.05, 0.85, 9.000, 0.950, 0.150, 0.900], + [3.500, 0.150, 0.45, 0.20, 4.500, 0.500, 0.800, 0.250], + [5.500, 0.050, 0.55, 0.10, 6.500, 0.650, 0.950, 0.400], + [4.000, 0.200, 0.40, 0.30, 5.000, 0.450, 0.700, 0.350], + ] +) + + +def _rebase(value, root: Path): + """Return the manifest with every input path anchored on the local copy. + + :param value: A manifest document, one of its sections or one of its values. + :param root: The directory the extracted datasets live in. + :type root: pathlib.Path + """ + if isinstance(value, dict): + return {key: _rebase(item, root) for key, item in value.items()} + if isinstance(value, str) and PATH_SEGMENT in value: + return str(root / value[value.index(PATH_SEGMENT) :]) + return value + + +def _input_paths(value, root: Path): + """Yield every path of a rebased document that lies under ``root``.""" + if isinstance(value, dict): + for item in value.values(): + yield from _input_paths(item, root) + elif isinstance(value, str) and value.startswith(str(root)): + yield value + + +def _read_evaluations(path: Path) -> list[dict[str, str]]: + """Return the rows of an ``evaluations.csv`` table, header first. + + A local copy of the reader of ``tests/unit/calibration/test_runner.py``: + importing one test module from another would run the module-level + ``importorskip`` of the other one during the collection of this one, for + three lines of CSV reading. + """ + with path.open(encoding="utf-8", newline="") as table: + return list(csv.DictReader(table)) + + +def _fixed_ldd(document: dict, destination: Path) -> str: + """Generate the drainage network of the basin once and return its path. + + The network is derived from the digital elevation model of the basin with + the parameters the model itself uses when no network is configured, and is + written so that every run of the calibration routes the runoff along the + very same network. + """ + import pcraster as pcr + + pcr.setclone(document["RASTERS"]["clone"]) + dem = pcr.readmap(document["RASTERS"]["dem"]) + pcr.report(pcr.lddcreate(dem, 1e31, 1e31, 1e31, 1e31), str(destination)) + return str(destination) + + +def test_the_calibration_recovers_the_parameters_of_the_ipojuca_basin( + dataset_dir, tmp_path, monkeypatch +): + manifest = dataset_dir / f"{BASIN}.json" + if not manifest.is_file(): + pytest.skip(f"{manifest} is not there") + + # The inputs are read where they were extracted, so the persistent auxiliary + # metadata of GDAL would annotate every raster of the dataset with a sidecar + # file of its statistics, next to the raster itself. The test writes nothing + # outside its own temporary directory; the workers inherit the variable. + monkeypatch.setenv("GDAL_PAM_ENABLED", "NO") + + assert [tuple(row) for row in INIT if not is_admissible(row)] == [], ( + "every member of the initial population must be admissible, or the search " + "would spend its generation on candidates it never runs" + ) + + document = _rebase(json.loads(manifest.read_text(encoding="utf-8")), dataset_dir) + missing = [path for path in _input_paths(document, dataset_dir) if not Path(path).exists()] + assert missing == [], "the manifest was rebased onto the local copy of the datasets" + + output_dir = tmp_path / "output" + output_dir.mkdir() + document["DIRECTORIES"]["output"] = str(output_dir) + document["SIM_TIME"] = {"start": "01/01/2000", "end": "01/06/2000"} + document["RASTERS"]["ldd"] = _fixed_ldd(document, tmp_path / "ldd.map") + + config_file = tmp_path / f"{BASIN}.json" + config_file.write_text(json.dumps(document, indent=2), encoding="utf-8") + parameters = ModelConfiguration(config_file, validate_input=False).calibration_parameters + expected = parameters.model_dump() + + written = Model.from_file(config_file).run().time_series["arn"][0] + observed = tmp_path / "observed.csv" + shutil.copyfile(written, observed) + assert len(observed.read_text(encoding="utf-8").splitlines()) == STEPS + 1, ( + "the header and one row per monthly step of the window" + ) + + temp_dir = tmp_path / "tmp" + result = calibrate( + config_file, + observed, + tmp_path / "run", + CalibrationSettings( + init=INIT, + maxiter=1, + popsize=5, + workers=2, + seed=1, + temp_dir=str(temp_dir), + ), + ) + + # The observed series is what the configuration itself wrote, so the global + # optimum of the search is that configuration, with an efficiency of 1. + assert result.best_objective == 0.0 + assert result.best_nse == pytest.approx(1.0, abs=1e-12) + assert result.best_parameters == pytest.approx(expected, abs=1e-12) + assert set(result.best_parameters) == set(CALIBRATION_PARAMETERS) + + rows = _read_evaluations(result.evaluations_csv) + assert list(rows[0]) == list(EVALUATION_COLUMNS) + assert len(rows) == result.evaluations + assert [row["error"] for row in rows if row["error"]] == [] + + summary = json.loads(result.result_json.read_text(encoding="utf-8")) + assert summary["nfev"] == result.evaluations + assert summary["population_size"] == len(INIT) + + assert result.calibrated_config.name == f"{BASIN}-calibrated.json" + calibrated = ModelConfiguration(result.calibrated_config, validate_input=False) + assert calibrated.calibration_parameters.model_dump() == pytest.approx( + result.best_parameters, abs=1e-12 + ) + + assert multiprocessing.active_children() == [] + assert list(temp_dir.iterdir()) == [] diff --git a/tests/unit/calibration/__init__.py b/tests/unit/calibration/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/unit/calibration/test_objective.py b/tests/unit/calibration/test_objective.py new file mode 100644 index 00000000..92ebc5f2 --- /dev/null +++ b/tests/unit/calibration/test_objective.py @@ -0,0 +1,695 @@ +"""The calibration objective: masks, Nash-Sutcliffe efficiency, tables and alignment.""" + +import dataclasses + +import numpy as np +import pytest + +from rubem.calibration.objective import ( + INADMISSIBLE_OBJECTIVE, + MISSING_VALUES, + PCRASTER_MISSING, + Series, + StationMetrics, + evaluate_series, + nash_sutcliffe, + objective, + read_series, + station_metrics, + valid_mask, +) + +# The PCRaster missing value as it comes back from a Float32 raster or table: +# 1e31 has no exact Float32 representation, so what is read is slightly below it. +FLOAT32_MISSING = float(np.float32(1e31)) + +# Every gap marker, used wherever a test drops a pair through one of them. +GAPS = (float("nan"), -9999.0, -1.0, 1e31, FLOAT32_MISSING) + + +def nse_of(metrics): + """The efficiency of every station of a mapping of :class:`StationMetrics`.""" + return {station: measured.nse for station, measured in metrics.items()} + + +def series(steps, **stations): + """A :class:`Series` built from plain lists, as the readers produce it.""" + return Series( + steps=np.asarray(steps, dtype=np.int64), + stations={ + station: np.asarray(values, dtype=np.float64) for station, values in stations.items() + }, + ) + + +class TestValidMask: + @pytest.mark.unit + def test_ordinary_values_are_valid(self): + assert valid_mask(np.array([0.0, 1.5, 3.25, 1e29])).tolist() == [True] * 4 + + @pytest.mark.unit + @pytest.mark.parametrize("negative", [-1.0, -999.0, -0.5, -9999.0]) + def test_a_negative_value_is_a_gap(self, negative): + # Every quantity the two series carry is a flux or a storage: a + # negative entry is a marker of a gap, whichever one a gauge writes. + assert valid_mask(np.array([negative])).tolist() == [False] + + @pytest.mark.unit + def test_zero_is_a_value_and_not_a_gap(self): + # A station may record no flow at all; only the sign below zero is a gap. + assert valid_mask(np.array([0.0, -0.0])).tolist() == [True, True] + + @pytest.mark.unit + @pytest.mark.parametrize("gap", GAPS) + def test_every_gap_marker_is_masked(self, gap): + assert valid_mask(np.array([1.0, gap, 2.0])).tolist() == [True, False, True] + + @pytest.mark.unit + def test_the_documented_markers_are_the_ones_the_mask_rejects(self): + assert MISSING_VALUES == (-9999.0,) + assert PCRASTER_MISSING == 1e31 + assert not valid_mask(np.array([PCRASTER_MISSING])).any() + + @pytest.mark.unit + def test_infinities_are_masked(self): + assert valid_mask(np.array([np.inf, -np.inf])).tolist() == [False, False] + + @pytest.mark.unit + def test_the_threshold_itself_is_the_first_value_rejected(self): + # The mask keeps everything strictly below 1e30 and drops 1e30 itself, + # which is what makes the Float32 round-trip of 1e31 a gap as well. + below = float(np.nextafter(1e30, 0.0)) + + assert valid_mask(np.array([below, 1e30])).tolist() == [True, False] + assert below < 1e30 < FLOAT32_MISSING + + +class TestNashSutcliffe: + @pytest.mark.unit + def test_an_exact_simulation_gives_one(self): + # residuals 0 / variability 2 -> NSE = 1 - 0 = 1.0 + observed = np.array([1.0, 2.0, 3.0]) + + assert nash_sutcliffe(observed.copy(), observed) == 1.0 + + @pytest.mark.unit + def test_the_mean_of_the_observations_gives_zero(self): + # observed [1, 2, 3], mean 2; simulated [2, 2, 2] + # residuals (1 + 0 + 1) = 2, variability (1 + 0 + 1) = 2 -> NSE = 1 - 2/2 = 0.0 + assert nash_sutcliffe(np.array([2.0, 2.0, 2.0]), np.array([1.0, 2.0, 3.0])) == 0.0 + + @pytest.mark.unit + def test_a_reversed_simulation_gives_a_negative_efficiency(self): + # observed [1, 2, 3], mean 2; simulated [3, 2, 1] + # residuals (4 + 0 + 4) = 8, variability 2 -> NSE = 1 - 8/2 = -3.0 + assert nash_sutcliffe(np.array([3.0, 2.0, 1.0]), np.array([1.0, 2.0, 3.0])) == -3.0 + + @pytest.mark.unit + @pytest.mark.parametrize("gap", GAPS) + def test_a_gap_in_the_observations_drops_the_pair(self, gap): + # The fourth pair is dropped, the first three are exact -> NSE = 1.0 + simulated = np.array([1.0, 2.0, 3.0, 99.0]) + observed = np.array([1.0, 2.0, 3.0, gap]) + + assert nash_sutcliffe(simulated, observed) == 1.0 + + @pytest.mark.unit + @pytest.mark.parametrize("gap", GAPS) + def test_a_gap_in_the_simulation_drops_the_pair(self, gap): + # Dropping the fourth pair leaves observed [1, 2, 3] with mean 2 and + # simulated [3, 2, 1]: residuals 8, variability 2 -> NSE = -3.0 + simulated = np.array([3.0, 2.0, 1.0, gap]) + observed = np.array([1.0, 2.0, 3.0, 4.0]) + + assert nash_sutcliffe(simulated, observed) == -3.0 + + @pytest.mark.unit + def test_a_gap_is_dropped_from_the_mean_of_the_observations(self): + # The third pair is dropped by the simulation, so the mean is that of + # [1, 3] = 2 and not that of [1, 3, 8]: variability (1 + 1) = 2, + # residuals (1 + 1) = 2 -> NSE = 1 - 2/2 = 0.0 + simulated = np.array([2.0, 2.0, float("nan")]) + observed = np.array([1.0, 3.0, 8.0]) + + assert nash_sutcliffe(simulated, observed) == 0.0 + + @pytest.mark.unit + def test_a_single_valid_pair_has_no_efficiency(self): + assert nash_sutcliffe(np.array([1.0, 2.0]), np.array([1.0, float("nan")])) is None + + @pytest.mark.unit + def test_a_station_without_one_valid_pair_has_no_efficiency(self): + # No pair survives the mask at all: the count guard must answer before + # anything is read out of the empty selection. + assert nash_sutcliffe(np.array([1.0, -9999.0]), np.array([np.nan, 2.0])) is None + assert nash_sutcliffe(np.array([np.nan, np.nan]), np.array([np.nan, np.nan])) is None + + @pytest.mark.unit + def test_two_valid_pairs_are_enough(self): + # observed [1, 3], mean 2, variability 2; simulated [1, 3] -> NSE = 1.0 + simulated = np.array([1.0, 3.0, 5.0]) + observed = np.array([1.0, 3.0, -9999.0]) + + assert nash_sutcliffe(simulated, observed) == 1.0 + + @pytest.mark.unit + def test_observations_without_variance_have_no_efficiency(self): + # The sum of the squared deviations of [0.1, 0.1, 0.1] is not exactly + # zero in floating point, so constancy is what the guard must test. + assert nash_sutcliffe(np.array([0.1, 0.2, 0.3]), np.array([0.1, 0.1, 0.1])) is None + + @pytest.mark.unit + def test_only_the_valid_observations_decide_the_variance(self): + # The valid observations are [5, 5]: constant, so there is no efficiency + # although the series as written also carries a 7. + assert nash_sutcliffe(np.array([5.0, 5.0, np.nan]), np.array([5.0, 5.0, 7.0])) is None + + @pytest.mark.unit + def test_series_of_different_lengths_are_refused(self): + with pytest.raises(ValueError, match="same shape"): + nash_sutcliffe(np.array([1.0, 2.0]), np.array([1.0, 2.0, 3.0])) + + +class TestObjective: + @pytest.mark.unit + @pytest.mark.parametrize( + ("nse", "expected"), + [ + # 1000 * (100 * (1 - 1))^2 = 0 + (1.0, 0.0), + # 1000 * (100 * 1)^2 = 1000 * 10000 = 1e7 + (0.0, 1.0e7), + # 1000 * (100 * 2)^2 = 1000 * 40000 = 4e7 + (-1.0, 4.0e7), + ], + ) + def test_the_objective_of_a_known_efficiency(self, nse, expected): + assert objective(nse) == expected + + @pytest.mark.unit + def test_the_inadmissible_objective_is_far_above_any_run(self): + # An NSE of -100 is an absurd simulation and still gives ~1.02e11. + assert objective(-100.0) < INADMISSIBLE_OBJECTIVE / 1e10 + assert np.isfinite(INADMISSIBLE_OBJECTIVE) + + +class TestReadSeries: + @pytest.mark.unit + def test_a_handwritten_table_in_the_model_csv_layout(self, tmp_path): + path = tmp_path / "tss_arn.csv" + path.write_text("0;1;2\n1;1.5;2.5\n2;3.5;4.5\n3;5.5;6.5\n", encoding="utf8") + + table = read_series(path) + + assert table.steps.tolist() == [1, 2, 3] + assert sorted(table.stations) == ["1", "2"] + assert table.stations["1"].tolist() == [1.5, 3.5, 5.5] + assert table.stations["2"].tolist() == [2.5, 4.5, 6.5] + + @pytest.mark.unit + def test_gaps_of_a_csv_table_are_kept_as_written(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;A;B\n1;-9999;2.0\n2;1e31;\n", encoding="utf8") + + table = read_series(path) + + assert table.stations["A"].tolist() == [-9999.0, 1e31] + assert table.stations["B"][0] == 2.0 + assert np.isnan(table.stations["B"][1]), "an empty cell is a gap" + assert valid_mask(table.stations["A"]).tolist() == [False, False] + + @pytest.mark.unit + def test_the_ids_of_a_real_run_are_the_ids_of_the_sample_raster(self, tmp_path): + """The table the model writes must be readable exactly as it lands.""" + from tests.unit.core.test_core import run_model + + run_model(str(tmp_path)) + + table = read_series(tmp_path / "out" / "tss_arn.csv") + + # The synthetic dataset marks two sample locations, with the ids 1 and 2. + assert sorted(table.stations) == ["1", "2"] + assert table.steps.tolist() == [1, 2] + for values in table.stations.values(): + assert values.shape == (2,) + assert valid_mask(values).all() + + @pytest.mark.unit + def test_the_windows_line_endings_the_model_writes(self, tmp_path): + # rubem.file._file_conversions.tss2csv writes through csv.writer, whose + # line terminator is CRLF; the reader must not take the carriage return + # for part of the last value of a row. + path = tmp_path / "tss_arn.csv" + path.write_bytes(b"0;1;2\r\n1;1.5;2.5\r\n2;3.5;4.5\r\n") + + table = read_series(path) + + assert table.steps.tolist() == [1, 2] + assert table.stations["2"].tolist() == [2.5, 4.5] + + @pytest.mark.unit + def test_the_exponent_spellings_of_the_pcraster_missing_value(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;A\n1;1e+31\n2;1E31\n3;2.0\n", encoding="utf8") + + table = read_series(path) + + assert valid_mask(table.stations["A"]).tolist() == [False, False, True] + + @pytest.mark.unit + def test_a_handwritten_pcraster_time_series(self, tmp_path): + path = tmp_path / "observed.tss" + path.write_text( + "Observed streamflow\n3\ntimestep\n11\n22\n1 1.5 2.5\n2 3.5 4.5\n", + encoding="utf8", + ) + + table = read_series(path) + + assert table.steps.tolist() == [1, 2] + assert sorted(table.stations) == ["11", "22"] + assert table.stations["11"].tolist() == [1.5, 3.5] + assert table.stations["22"].tolist() == [2.5, 4.5] + + @pytest.mark.unit + def test_a_pcraster_header_whose_count_carries_trailing_spaces(self, tmp_path): + # The number of columns is read from a line of its own, which a writer + # may pad; the title above it carries spaces of its own and is not read. + path = tmp_path / "observed.tss" + path.write_text( + "timeseries scalar\n3 \ntimestep\n11\n22\n1 1.5 2.5\n2 3.5 4.5\n", + encoding="utf8", + ) + + table = read_series(path) + + assert sorted(table.stations) == ["11", "22"] + assert table.stations["22"].tolist() == [2.5, 4.5] + + @pytest.mark.unit + def test_the_windows_line_endings_of_a_pcraster_time_series(self, tmp_path): + # The carriage return must not be taken for part of the number of + # columns, of a station id or of the last value of a row. + path = tmp_path / "observed.tss" + path.write_bytes( + b"timeseries scalar\r\n3\r\ntimestep\r\n11\r\n22\r\n1 1.5 2.5\r\n2 3.5 4.5\r\n" + ) + + table = read_series(path) + + assert sorted(table.stations) == ["11", "22"] + assert table.steps.tolist() == [1, 2] + assert table.stations["22"].tolist() == [2.5, 4.5] + + @pytest.mark.unit + def test_a_pcraster_time_series_without_the_header_is_refused(self, tmp_path): + # Numbering the columns 1..N would label them with a guess: the columns + # of a time series are stations, and nothing here says which. + path = tmp_path / "observed.tss" + path.write_text("1 1.5 2.5\n2 3.5 4.5\n", encoding="utf8") + + with pytest.raises(ValueError, match="PCRaster time series files carry") as failure: + read_series(path) + + message = str(failure.value) + assert "has no header" in message + assert "0;;..." in message, "the message names both accepted layouts" + assert "title line" in message + + @pytest.mark.unit + @pytest.mark.parametrize( + "content", + [ + "Observed streamflow\n1\ntimestep\n1 1.5\n", + "Observed streamflow\n3\ntimestep\n11\n", + "1 1.5 2.5\n", + ], + ids=["one column", "fewer names than the count", "a single row"], + ) + def test_a_truncated_pcraster_header_is_refused(self, tmp_path, content): + path = tmp_path / "observed.tss" + path.write_text(content, encoding="utf8") + + with pytest.raises(ValueError, match="PCRaster time series files carry"): + read_series(path) + + @pytest.mark.unit + def test_a_csv_table_without_the_header_is_refused(self, tmp_path): + # The first record is already a time step, so the ids would be read + # from it and the first step would be lost from the series. + path = tmp_path / "observed.csv" + path.write_text("1;16.33;22.94\n2;14.86;19.95\n", encoding="utf8") + + with pytest.raises(ValueError, match="has no header") as failure: + read_series(path) + + message = str(failure.value) + assert "'1'" in message, "the message quotes the line it refused" + assert "0;;..." in message, "the message names both accepted layouts" + assert "PCRaster time series" in message + + @pytest.mark.unit + def test_a_byte_order_mark_is_not_part_of_the_step_column_label(self, tmp_path): + # An editor that writes UTF-8 with a byte order mark puts it at the + # start of the first line, where it would otherwise become part of the + # first header cell. + path = tmp_path / "observed.csv" + path.write_bytes("0;A;B\n1;1.5;2.5\n2;3.5;4.5\n".encode("utf8")) + + table = read_series(path) + + assert sorted(table.stations) == ["A", "B"] + assert table.steps.tolist() == [1, 2] + assert table.stations["A"].tolist() == [1.5, 3.5] + + @pytest.mark.unit + def test_a_byte_order_mark_does_not_hide_a_missing_csv_header(self, tmp_path): + # The mark would turn the first time step into a word, and a word is + # the label of a step column: the refusal must not depend on it. + path = tmp_path / "observed.csv" + path.write_bytes("1;16.33;22.94\n2;14.86;19.95\n".encode("utf8")) + + with pytest.raises(ValueError, match="has no header") as failure: + read_series(path) + + assert "'1'" in str(failure.value), "the mark is not part of the line it quotes" + + @pytest.mark.unit + def test_a_byte_order_mark_before_the_title_of_a_time_series(self, tmp_path): + path = tmp_path / "observed.tss" + path.write_bytes( + "timeseries scalar\n3\ntimestep\n11\n22\n1 1.5 2.5\n2 3.5 4.5\n".encode("utf8") + ) + + table = read_series(path) + + assert sorted(table.stations) == ["11", "22"] + assert table.steps.tolist() == [1, 2] + + @pytest.mark.unit + @pytest.mark.parametrize("label", ["0", "timestep", "step", ""]) + def test_the_accepted_spellings_of_the_step_column_label(self, tmp_path, label): + path = tmp_path / "observed.csv" + path.write_text(f"{label};A\n1;1.0\n2;2.0\n", encoding="utf8") + + table = read_series(path) + + assert table.steps.tolist() == [1, 2] + assert table.stations["A"].tolist() == [1.0, 2.0] + + @pytest.mark.unit + def test_a_row_that_does_not_start_with_a_time_step_names_the_file(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;A\n2000-01-01;1.0\n2000-01-02;2.0\n", encoding="utf8") + + with pytest.raises(ValueError, match="Row 2 .* does not start with a time step"): + read_series(path) + + @pytest.mark.unit + def test_a_repeated_station_id_is_refused(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;1;1\n1;1.0;2.0\n", encoding="utf8") + + with pytest.raises(ValueError, match="repeated station ids: 1"): + read_series(path) + + @pytest.mark.unit + def test_a_row_that_does_not_match_the_header_is_refused(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;1;2\n1;1.0;2.0\n2;3.0\n", encoding="utf8") + + with pytest.raises(ValueError, match="Row 3"): + read_series(path) + + @pytest.mark.unit + def test_an_empty_table_is_refused(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("", encoding="utf8") + + with pytest.raises(ValueError, match="empty"): + read_series(path) + + +class TestEvaluateSeries: + @pytest.mark.unit + def test_the_mean_over_the_stations(self): + # Station A: exact -> 1.0. Station B: observed [1, 2, 3] with the + # simulation reversed -> -3.0. Mean = (1.0 - 3.0) / 2 = -1.0 + simulated = series([1, 2, 3], A=[4.0, 5.0, 6.0], B=[3.0, 2.0, 1.0]) + observed = series([1, 2, 3], A=[4.0, 5.0, 6.0], B=[1.0, 2.0, 3.0]) + + mean, per_station = evaluate_series(simulated, observed) + + assert nse_of(per_station) == {"A": 1.0, "B": -3.0} + assert mean == -1.0 + + @pytest.mark.unit + def test_only_the_shared_steps_are_compared(self): + # The shared steps are 3 and 4: observed [1, 3] with mean 2 and + # variability 2, simulated [1, 3] -> NSE = 1.0 + simulated = series([1, 2, 3, 4], A=[99.0, 99.0, 1.0, 3.0]) + observed = series([3, 4, 5], A=[1.0, 3.0, 99.0]) + + mean, per_station = evaluate_series(simulated, observed) + + assert nse_of(per_station) == {"A": 1.0} + assert per_station["A"].pairs == 2, "only the shared steps are a pair" + assert mean == 1.0 + + @pytest.mark.unit + def test_the_spin_up_steps_are_excluded(self): + # Without the spin-up the first step would make the simulation exact; + # with two spin-up steps only the steps 3 and 4 are compared: + # observed [1, 3], mean 2, variability 2; simulated [3, 1] + # -> residuals (4 + 4) = 8 -> NSE = 1 - 8/2 = -3.0 + simulated = series([1, 2, 3, 4], A=[0.0, 0.0, 3.0, 1.0]) + observed = series([1, 2, 3, 4], A=[0.0, 0.0, 1.0, 3.0]) + + mean, per_station = evaluate_series(simulated, observed, spinup_steps=2) + + assert nse_of(per_station) == {"A": -3.0} + assert mean == -3.0 + + @pytest.mark.unit + def test_a_station_present_on_one_side_only_is_ignored(self): + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[1.0, 2.0, 3.0]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0], C=[1.0, 2.0, 3.0]) + + mean, per_station = evaluate_series(simulated, observed) + + assert nse_of(per_station) == {"A": 1.0} + assert mean == 1.0 + + @pytest.mark.unit + def test_a_station_without_variance_is_reported_and_left_out_of_the_mean(self): + # Station A is exact -> 1.0; station B has constant observations -> None. + # The mean is that of the stations with a value: 1.0 + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[0.1, 0.2, 0.3]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[0.1, 0.1, 0.1]) + + mean, per_station = evaluate_series(simulated, observed) + + assert nse_of(per_station) == {"A": 1.0, "B": None} + assert per_station["B"].pairs == 3, "the station is measured even without an efficiency" + assert mean == 1.0 + + @pytest.mark.unit + def test_no_shared_station_names_both_sides(self): + simulated = series([1, 2], A=[1.0, 2.0]) + observed = series([1, 2], B=[1.0, 2.0]) + + with pytest.raises(ValueError, match="no station in common") as error: + evaluate_series(simulated, observed) + + assert "A" in str(error.value) and "B" in str(error.value) + + @pytest.mark.unit + def test_no_shared_step_names_the_steps(self): + simulated = series([1, 2], A=[1.0, 2.0]) + observed = series([7, 8], A=[1.0, 2.0]) + + with pytest.raises(ValueError, match="no time step in common") as error: + evaluate_series(simulated, observed) + + assert "[1, 2]" in str(error.value) and "[7, 8]" in str(error.value) + + @pytest.mark.unit + def test_a_spin_up_that_eats_every_step_is_an_error(self): + simulated = series([1, 2], A=[1.0, 2.0]) + observed = series([1, 2], A=[1.0, 2.0]) + + with pytest.raises(ValueError, match="no time step in common"): + evaluate_series(simulated, observed, spinup_steps=2) + + @pytest.mark.unit + def test_no_station_with_an_efficiency_is_an_error(self): + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[1.0, 2.0, 3.0]) + observed = series([1, 2, 3], A=[5.0, 5.0, 5.0], B=[1.0, -9999.0, np.nan]) + + with pytest.raises(ValueError, match="No station has a Nash-Sutcliffe") as error: + evaluate_series(simulated, observed) + + assert "A, B" in str(error.value) + + +class TestStationMetrics: + @pytest.mark.unit + def test_every_statistic_of_a_hand_computed_pair_of_series(self): + # observed [1, 2, 6], simulated [2, 2, 5]. + # mean_observed = (1 + 2 + 6) / 3 = 3 + # mean_simulated = (2 + 2 + 5) / 3 = 3 + # deviations: observed [-2, -1, 3], simulated [-1, -1, 2] + # std_observed = sqrt((4 + 1 + 9) / 2) = sqrt(7) + # std_simulated = sqrt((1 + 1 + 4) / 2) = sqrt(3) + # r = (2 + 1 + 6) / sqrt(14 * 6) = 9 / sqrt(84) + # rmse = sqrt(((2-1)^2 + 0 + (5-6)^2) / 3) = sqrt(2/3) + # nse = 1 - (1 + 0 + 1) / 14 = 1 - 2/14 = 6/7 + measured = station_metrics(np.array([2.0, 2.0, 5.0]), np.array([1.0, 2.0, 6.0])) + + assert measured.pairs == 3 + assert measured.mean_observed == pytest.approx(3.0) + assert measured.mean_simulated == pytest.approx(3.0) + assert measured.std_observed == pytest.approx(np.sqrt(7.0)) + assert measured.std_simulated == pytest.approx(np.sqrt(3.0)) + assert measured.r == pytest.approx(9.0 / np.sqrt(84.0)) + assert measured.rmse == pytest.approx(np.sqrt(2.0 / 3.0)) + assert measured.nse == pytest.approx(6.0 / 7.0) + + @pytest.mark.unit + def test_it_is_a_frozen_record_of_eight_fields(self): + measured = station_metrics(np.array([1.0, 2.0]), np.array([1.0, 2.0])) + + assert isinstance(measured, StationMetrics) + assert list(dataclasses.asdict(measured)) == [ + "pairs", + "mean_observed", + "std_observed", + "mean_simulated", + "std_simulated", + "r", + "rmse", + "nse", + ] + with pytest.raises(dataclasses.FrozenInstanceError): + measured.nse = 0.0 + + @pytest.mark.unit + def test_only_the_pairs_valid_in_both_series_are_measured(self): + # The third step is a gap in the simulation and the fourth a negative + # observation: the sample is observed [1, 3] and simulated [1, 3]. + # mean 2 on both sides, std sqrt(((1)^2 + (1)^2) / 1) = sqrt(2) + measured = station_metrics( + np.array([1.0, 3.0, np.nan, 9.0]), np.array([1.0, 3.0, 8.0, -1.0]) + ) + + assert measured.pairs == 2 + assert measured.mean_observed == pytest.approx(2.0) + assert measured.std_observed == pytest.approx(np.sqrt(2.0)) + assert measured.rmse == 0.0 + assert measured.nse == 1.0 + + @pytest.mark.unit + def test_a_single_pair_has_a_mean_and_an_error_but_no_spread(self): + # One pair: the mean and the root mean squared error are defined, a + # standard deviation with one degree of freedom and an efficiency are not. + measured = station_metrics(np.array([3.0, np.nan]), np.array([1.0, 2.0])) + + assert measured.pairs == 1 + assert measured.mean_observed == 1.0 + assert measured.mean_simulated == 3.0 + assert measured.rmse == 2.0 + assert measured.std_observed is None + assert measured.std_simulated is None + assert measured.r is None + assert measured.nse is None + + @pytest.mark.unit + def test_a_station_without_one_valid_pair_has_nothing(self): + measured = station_metrics(np.array([1.0, -9999.0]), np.array([np.nan, 2.0])) + + assert measured.pairs == 0 + assert dataclasses.astuple(measured) == (0, None, None, None, None, None, None, None) + + @pytest.mark.unit + def test_a_constant_series_has_no_correlation(self): + # The simulation never moves, so its own spread is zero and the + # correlation has no denominator; the efficiency is still defined, + # since the observations vary. + measured = station_metrics(np.array([2.0, 2.0, 2.0]), np.array([1.0, 2.0, 3.0])) + + assert measured.std_simulated == 0.0 + assert measured.r is None + assert measured.nse == 0.0 + + @pytest.mark.unit + def test_constant_observations_leave_neither_a_correlation_nor_an_efficiency(self): + measured = station_metrics(np.array([0.1, 0.2, 0.3]), np.array([0.1, 0.1, 0.1])) + + # The spread of a constant series is only approximately zero in + # floating point, which is why the guards test the values themselves. + assert measured.std_observed == pytest.approx(0.0, abs=1e-15) + assert measured.r is None + assert measured.nse is None + + @pytest.mark.unit + def test_series_of_different_lengths_are_refused(self): + with pytest.raises(ValueError, match="same shape"): + station_metrics(np.array([1.0, 2.0]), np.array([1.0, 2.0, 3.0])) + + +class TestStationSelection: + @pytest.mark.unit + def test_the_mean_is_taken_over_the_selection_and_every_station_is_measured(self): + # Station A is exact -> 1.0, station B reversed -> -3.0, station C + # exact -> 1.0. The objective averages A and C: (1.0 + 1.0) / 2 = 1.0, + # and B is measured all the same. + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[3.0, 2.0, 1.0], C=[4.0, 5.0, 6.0]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[1.0, 2.0, 3.0], C=[4.0, 5.0, 6.0]) + + mean, per_station = evaluate_series(simulated, observed, stations=("A", "C")) + + assert mean == 1.0 + assert nse_of(per_station) == {"A": 1.0, "B": -3.0, "C": 1.0} + + @pytest.mark.unit + def test_a_selected_station_without_an_efficiency_is_left_out_of_the_mean(self): + # B has constant observations and no efficiency, so the mean of the + # selection is that of A alone. + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[0.1, 0.2, 0.3]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[0.1, 0.1, 0.1]) + + mean, _ = evaluate_series(simulated, observed, stations=("A", "B")) + + assert mean == 1.0 + + @pytest.mark.unit + def test_a_station_of_the_selection_the_series_do_not_share_is_ignored(self): + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[3.0, 2.0, 1.0]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[1.0, 2.0, 3.0]) + + mean, per_station = evaluate_series(simulated, observed, stations=("A", "Z")) + + assert mean == 1.0 + assert sorted(per_station) == ["A", "B"] + + @pytest.mark.unit + def test_a_selection_that_names_no_shared_station_names_the_ids(self): + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0]) + + with pytest.raises(ValueError, match="none of the stations") as failure: + evaluate_series(simulated, observed, stations=("Y", "Z")) + + message = str(failure.value) + assert "Y, Z" in message + assert "A" in message + + @pytest.mark.unit + def test_no_selected_station_with_an_efficiency_names_the_selection(self): + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[1.0, 2.0, 3.0]) + observed = series([1, 2, 3], A=[5.0, 5.0, 5.0], B=[1.0, 2.0, 3.0]) + + with pytest.raises(ValueError, match="No station has a Nash-Sutcliffe") as failure: + evaluate_series(simulated, observed, stations=("A",)) + + assert "(A)" in str(failure.value), "only the selection is named" diff --git a/tests/unit/calibration/test_parameters.py b/tests/unit/calibration/test_parameters.py new file mode 100644 index 00000000..a2c209a7 --- /dev/null +++ b/tests/unit/calibration/test_parameters.py @@ -0,0 +1,455 @@ +"""The decision vector of the calibration: bounds, conversions and admissibility.""" + +import numpy as np +import pytest + +from rubem.calibration.parameters import ( + CALIBRATION_PARAMETERS, + DERIVED_PARAMETER, + FREE_PARAMETERS, + bounds, + decision_space, + is_admissible, + parameters_to_vector, + vector_to_parameters, + weights_constraint, +) +from rubem.configuration.app_settings import AppSettings +from rubem.configuration.calibration_parameters import CalibrationParameters +from rubem.configuration.model_configuration import ModelConfiguration +from tests.helpers.synthetic import write_synthetic_dataset + +# The ranges of rubem/appsettings.json, written out so that a change of the +# settings file, or of the order of the decision vector, has to be made here +# too instead of passing unnoticed through a comparison with the same source. +SETTINGS_BOUNDS = { + "alpha": (0.01, 10.0), + "beta": (0.01, 1.0), + "w_1": (0.0, 1.0), + "w_2": (0.0, 1.0), + "w_3": (0.0, 1.0), + "rcd": (1.0, 10.0), + "f": (0.01, 1.0), + "alpha_gw": (0.01, 1.0), + "x": (0.0, 1.0), +} + + +def admissible_vector(): + """A decision vector well inside every bound, with weights that add up below 1.""" + return np.array([4.5, 0.5, 0.333, 0.333, 5.0, 0.5, 0.5, 0.5], dtype=np.float64) + + +class TestBounds: + @pytest.mark.unit + def test_the_bounds_are_the_ranges_of_the_application_settings(self): + ranges = AppSettings.default().value_ranges.variables + + assert bounds() == [(ranges[name].min, ranges[name].max) for name in FREE_PARAMETERS] + + @pytest.mark.unit + def test_the_bounds_are_the_ones_written_in_the_settings_file(self): + assert bounds() == [SETTINGS_BOUNDS[name] for name in FREE_PARAMETERS] + assert [name for name in SETTINGS_BOUNDS if name != DERIVED_PARAMETER] == list( + FREE_PARAMETERS + ) + + @pytest.mark.unit + def test_the_free_parameters_do_not_include_the_derived_weight(self): + assert DERIVED_PARAMETER not in FREE_PARAMETERS + assert set(CALIBRATION_PARAMETERS) == set(FREE_PARAMETERS) | {DERIVED_PARAMETER} + assert len(bounds()) == len(FREE_PARAMETERS) == 8 + + +class TestConversions: + @pytest.mark.unit + def test_a_vector_becomes_the_nine_parameters_with_the_derived_weight(self): + parameters = vector_to_parameters(admissible_vector()) + + assert list(parameters) == list(CALIBRATION_PARAMETERS) + assert parameters["beta"] == 0.5 + assert parameters["w_3"] == 1.0 - (0.333 + 0.333) + assert parameters["w_1"] + parameters["w_2"] + parameters["w_3"] == pytest.approx(1.0) + + @pytest.mark.unit + def test_a_vector_of_the_wrong_size_is_rejected(self): + with pytest.raises(ValueError, match="must have 8 entries"): + vector_to_parameters([1.0, 2.0, 3.0]) + + @pytest.mark.unit + def test_the_configuration_parameters_round_trip_through_the_vector(self, tmp_path): + config = write_synthetic_dataset(str(tmp_path)) + configuration = ModelConfiguration(config, validate_input=False) + parameters = configuration.calibration_parameters.model_dump() + + vector = parameters_to_vector(parameters) + + assert vector_to_parameters(vector) == pytest.approx(parameters, abs=1e-12) + + @pytest.mark.unit + def test_the_legacy_spellings_of_the_parameters_are_accepted(self): + legacy = { + "alpha": 4.5, + "b": 0.5, + "w1": 0.333, + "w2": 0.333, + "w3": 0.334, + "rcd": 5.0, + "f": 0.5, + "alpha_gw": 0.5, + "x": 0.5, + } + + assert parameters_to_vector(legacy) == pytest.approx(admissible_vector()) + + @pytest.mark.unit + def test_a_missing_parameter_names_what_the_vector_needs(self): + with pytest.raises(KeyError, match="alpha_gw"): + parameters_to_vector({"alpha": 4.5, "beta": 0.5}) + + +class TestAdmissibility: + @pytest.mark.unit + def test_a_vector_inside_the_bounds_is_admissible(self): + assert is_admissible(admissible_vector()) + + @pytest.mark.unit + def test_a_vector_on_the_bounds_is_admissible(self): + vector = np.array([minimum for minimum, _ in bounds()], dtype=np.float64) + + assert is_admissible(vector) + + @pytest.mark.unit + @pytest.mark.parametrize("index", range(8)) + def test_a_value_above_its_bound_is_rejected(self, index): + vector = admissible_vector() + vector[index] = bounds()[index][1] + 1.0 + + assert not is_admissible(vector) + + @pytest.mark.unit + @pytest.mark.parametrize("index", range(8)) + def test_a_value_below_its_bound_is_rejected(self, index): + vector = admissible_vector() + vector[index] = bounds()[index][0] - 1.0 + + assert not is_admissible(vector) + + @pytest.mark.unit + def test_weights_that_add_up_above_one_are_rejected(self): + vector = admissible_vector() + vector[FREE_PARAMETERS.index("w_1")] = 0.6 + vector[FREE_PARAMETERS.index("w_2")] = 0.6 + + # Both weights are inside their own [0, 1] bound; what makes the + # candidate inadmissible is their sum, and the negative w_3 it derives. + assert vector_to_parameters(vector)["w_3"] < 0.0 + assert not is_admissible(vector) + + @pytest.mark.unit + def test_weights_that_add_up_to_exactly_one_are_admissible(self): + vector = admissible_vector() + vector[FREE_PARAMETERS.index("w_1")] = 0.4 + vector[FREE_PARAMETERS.index("w_2")] = 0.6 + + assert vector_to_parameters(vector)["w_3"] == 0.0 + assert is_admissible(vector) + + @pytest.mark.unit + @pytest.mark.parametrize("hundredths", range(101)) + def test_every_pair_of_weights_that_adds_up_to_one_reaches_the_model(self, hundredths): + """What the guard admits, the configuration model must accept. + + ``1 - w_1 - w_2`` rounds twice and lands one ulp below zero for 20 of + these 101 pairs, ``w_1 = 0.33, w_2 = 0.67`` among them, which is a + weight :class:`CalibrationParameters` refuses; deriving the weight with + a single subtraction gives exactly zero for all of them. + """ + vector = admissible_vector() + vector[FREE_PARAMETERS.index("w_1")] = hundredths / 100.0 + vector[FREE_PARAMETERS.index("w_2")] = (100 - hundredths) / 100.0 + parameters = vector_to_parameters(vector) + + assert parameters["w_3"] == 0.0 + assert is_admissible(vector) + # No exception: the candidate the guard admits is one the model runs. + CalibrationParameters(**parameters) + + @pytest.mark.unit + def test_the_guard_and_the_derivation_cannot_disagree(self): + """Every admissible candidate of a grid is one the model accepts.""" + refused = [] + for first in range(0, 101, 7): + for second in range(0, 101, 11): + vector = admissible_vector() + vector[FREE_PARAMETERS.index("w_1")] = first / 100.0 + vector[FREE_PARAMETERS.index("w_2")] = second / 100.0 + if not is_admissible(vector): + continue + try: + CalibrationParameters(**vector_to_parameters(vector)) + except ValueError: + refused.append((first / 100.0, second / 100.0)) + + assert not refused, f"admitted but refused by the model: {refused}" + + @pytest.mark.unit + def test_a_vector_of_the_wrong_size_or_with_a_nan_is_rejected(self): + vector = admissible_vector() + vector[0] = np.nan + + assert not is_admissible(vector) + assert not is_admissible(admissible_vector()[:-1]) + + +class TestWeightsConstraint: + """The only part of the module that needs the optional SciPy extra.""" + + @pytest.fixture(autouse=True) + def _scipy(self): + pytest.importorskip("scipy.optimize") + + @pytest.mark.unit + def test_the_constraint_is_the_sum_of_the_two_searched_weights(self): + constraint = weights_constraint() + + expected = np.zeros((1, len(FREE_PARAMETERS))) + expected[0, FREE_PARAMETERS.index("w_1")] = 1.0 + expected[0, FREE_PARAMETERS.index("w_2")] = 1.0 + assert np.array_equal(np.asarray(constraint.A), expected) + assert np.all(np.asarray(constraint.lb) == -np.inf) + assert np.all(np.asarray(constraint.ub) == 1.0) + + @pytest.mark.unit + def test_the_constraint_evaluates_to_the_sum_of_the_weights(self): + constraint = weights_constraint() + vector = admissible_vector() + + assert np.asarray(constraint.A) @ vector == pytest.approx([0.666]) + + +class TestDefaultDecisionSpace: + @pytest.mark.unit + def test_a_space_without_arguments_is_the_whole_search(self): + space = decision_space() + + assert space.free_names == FREE_PARAMETERS + assert space.fixed == {} + assert space.dimension == 8 + assert list(space.bounds) == bounds() + + @pytest.mark.unit + def test_the_module_level_functions_are_the_default_space(self): + space = decision_space() + vector = admissible_vector() + parameters = vector_to_parameters(vector) + + assert space.to_parameters(vector) == parameters + assert space.from_parameters(parameters) == pytest.approx(parameters_to_vector(parameters)) + assert space.is_admissible(vector) is is_admissible(vector) is True + assert np.array_equal( + np.asarray(space.weights_constraint().A), np.asarray(weights_constraint().A) + ) + + +class TestFixedParameters: + @pytest.mark.unit + def test_a_fixed_parameter_leaves_the_vector_and_stays_in_the_parameters(self): + space = decision_space(fixed={"x": 0.0}) + + assert space.free_names == tuple(name for name in FREE_PARAMETERS if name != "x") + assert space.dimension == 7 + assert space.fixed == {"x": 0.0} + assert list(space.bounds) == [ + SETTINGS_BOUNDS[name] for name in FREE_PARAMETERS if name != "x" + ] + + parameters = space.to_parameters(np.delete(admissible_vector(), FREE_PARAMETERS.index("x"))) + + assert list(parameters) == list(CALIBRATION_PARAMETERS) + assert parameters["x"] == 0.0 + assert parameters["alpha"] == 4.5 + assert parameters["w_3"] == pytest.approx(1.0 - 0.666) + + @pytest.mark.unit + def test_the_free_vector_of_a_space_with_fixed_parameters_round_trips(self): + space = decision_space(fixed={"x": 0.25, "alpha": 2.0}) + parameters = vector_to_parameters(admissible_vector()) + + vector = space.from_parameters(parameters) + + assert len(vector) == space.dimension == 6 + # The fixed coordinates are dropped on the way in and restored, at the + # value they were pinned to, on the way out. + assert space.to_parameters(vector) == pytest.approx( + {**parameters, "x": 0.25, "alpha": 2.0}, abs=1e-12 + ) + + @pytest.mark.unit + def test_the_legacy_spellings_name_the_parameter_that_is_fixed(self): + space = decision_space(fixed={"b": 0.25, "w1": 0.5}) + + assert space.fixed == {"beta": 0.25, "w_1": 0.5} + assert "beta" not in space.free_names + + @pytest.mark.unit + def test_a_vector_of_the_wrong_size_names_the_free_parameters(self): + space = decision_space(fixed={"x": 0.0}) + + with pytest.raises(ValueError, match="must have 7 entries"): + space.to_parameters(admissible_vector()) + + @pytest.mark.unit + def test_the_derived_weight_cannot_be_fixed(self): + with pytest.raises(ValueError, match="w_3"): + decision_space(fixed={"w_3": 0.3}) + with pytest.raises(ValueError, match="derived"): + decision_space(fixed={"w3": 0.3}) + + @pytest.mark.unit + def test_an_unknown_parameter_cannot_be_fixed(self): + with pytest.raises(ValueError, match="not a searched calibration parameter"): + decision_space(fixed={"gamma": 1.0}) + + @pytest.mark.unit + @pytest.mark.parametrize("value", [-0.5, 1.5, float("nan")]) + def test_a_fixed_value_outside_the_settings_range_is_refused(self, value): + with pytest.raises(ValueError, match=r"'x' cannot be fixed"): + decision_space(fixed={"x": value}) + + @pytest.mark.unit + def test_fixing_every_parameter_leaves_nothing_to_search(self): + every = {name: minimum for name, (minimum, _) in SETTINGS_BOUNDS.items()} + every.pop(DERIVED_PARAMETER) + + with pytest.raises(ValueError, match="nothing to look for"): + decision_space(fixed=every) + + +class TestFixedWeights: + @pytest.mark.unit + def test_a_weight_left_with_one_value_is_fixed_as_well(self): + """w_1 = 1 leaves w_2 no room: w_2 = 0 is fixed and leaves the vector. + + The dimension, the budget and the initial population then describe the + search that runs, and no zero-width coordinate reaches the optimizer. + """ + space = decision_space(fixed={"w_1": 1.0}) + + assert "w_2" not in space.free_names + assert space.fixed == {"w_1": 1.0, "w_2": 0.0} + assert space.dimension == 6 + assert space.weights_constraint() is None + parameters = space.to_parameters([4.5, 0.5, 5.0, 0.5, 0.5, 0.5]) + assert (parameters["w_1"], parameters["w_2"], parameters["w_3"]) == (1.0, 0.0, 0.0) + + @pytest.mark.unit + def test_fixing_one_weight_narrows_the_other_instead_of_constraining_the_pair(self): + space = decision_space(fixed={"w_1": 0.7}) + position = space.free_names.index("w_2") + + # w_3 = 1 - w_1 - w_2 may not fall below its own minimum, so w_2 is + # searched in [0, 0.3] and the pair needs no constraint of its own. + assert space.bounds[position] == (0.0, pytest.approx(0.3)) + assert space.weights_constraint() is None + + @pytest.mark.unit + def test_the_narrowed_bound_keeps_the_override_when_the_override_is_tighter(self): + space = decision_space(fixed={"w_2": 0.5}, bounds={"w_1": (0.1, 0.2)}) + position = space.free_names.index("w_1") + + assert space.bounds[position] == (0.1, 0.2) + + @pytest.mark.unit + def test_a_fixed_weight_that_leaves_the_other_no_value_is_refused(self): + with pytest.raises(ValueError, match="no value in its range"): + decision_space(fixed={"w_1": 0.9}, bounds={"w_2": (0.2, 0.5)}) + + @pytest.mark.unit + def test_fixing_both_weights_checks_the_derived_one(self): + space = decision_space(fixed={"w_1": 0.3, "w_2": 0.4}) + vector = np.array([4.5, 0.5, 5.0, 0.5, 0.5, 0.5], dtype=np.float64) + + assert space.free_names == ("alpha", "beta", "rcd", "f", "alpha_gw", "x") + assert space.weights_constraint() is None + assert space.to_parameters(vector)["w_3"] == pytest.approx(0.3) + assert space.is_admissible(vector) + + @pytest.mark.unit + def test_two_fixed_weights_that_add_up_above_one_are_refused(self): + with pytest.raises(ValueError, match="derives w_3"): + decision_space(fixed={"w_1": 0.6, "w_2": 0.6}) + + @pytest.mark.unit + def test_the_constraint_follows_the_positions_of_the_weights_in_the_vector(self): + pytest.importorskip("scipy.optimize") + space = decision_space(fixed={"alpha": 5.0}) + constraint = space.weights_constraint() + + expected = np.zeros((1, space.dimension)) + expected[0, space.free_names.index("w_1")] = 1.0 + expected[0, space.free_names.index("w_2")] = 1.0 + assert np.array_equal(np.asarray(constraint.A), expected) + assert np.all(np.asarray(constraint.ub) == 1.0) + + +class TestOverriddenBounds: + @pytest.mark.unit + def test_an_override_narrows_the_range_of_the_settings(self): + space = decision_space(bounds={"rcd": (2.0, 5.0)}) + position = space.free_names.index("rcd") + + assert space.bounds[position] == (2.0, 5.0) + assert [ + bound + for name, bound in zip(space.free_names, space.bounds, strict=True) + if name != "rcd" + ] == [SETTINGS_BOUNDS[name] for name in FREE_PARAMETERS if name != "rcd"] + + @pytest.mark.unit + def test_a_candidate_outside_the_narrowed_bound_is_not_admissible(self): + space = decision_space(bounds={"rcd": (2.0, 5.0)}) + vector = admissible_vector() + + assert space.is_admissible(vector) + vector[FREE_PARAMETERS.index("rcd")] = 8.0 + assert not space.is_admissible(vector) + # The same candidate is admissible in the default space: the refusal is + # the narrowed bound of this run and not a bound of the model. + assert is_admissible(vector) + + @pytest.mark.unit + @pytest.mark.parametrize( + "override", + [(0.5, 5.0), (2.0, 12.0), (5.0, 2.0), (3.0, 3.0), (float("nan"), 5.0)], + ids=["below", "above", "reversed", "empty", "not a number"], + ) + def test_a_bound_that_does_not_narrow_the_settings_range_is_refused(self, override): + with pytest.raises(ValueError, match=r"of 'rcd' is not a narrower range") as failure: + decision_space(bounds={"rcd": override}) + + # Both ranges are named, so the message says what was asked and what is + # available. + assert "(1.0, 10.0)" in str(failure.value) + + @pytest.mark.unit + @pytest.mark.parametrize( + "override", [3.0, (3.0,), (1.0, 2.0, 3.0)], ids=["a number", "one value", "three values"] + ) + def test_a_bound_that_is_not_a_pair_is_refused(self, override): + # A bound is a range: anything else is answered with the pair it needs + # and not with the error of the conversion it failed. + with pytest.raises(ValueError, match=r"bound of 'rcd' must be a \(minimum, maximum\) pair"): + decision_space(bounds={"rcd": override}) + + @pytest.mark.unit + def test_a_bound_of_the_derived_or_of_an_unknown_parameter_is_refused(self): + with pytest.raises(ValueError, match="derived"): + decision_space(bounds={"w_3": (0.1, 0.5)}) + with pytest.raises(ValueError, match="not a searched calibration parameter"): + decision_space(bounds={"gamma": (0.1, 0.5)}) + + @pytest.mark.unit + def test_a_bound_of_a_fixed_parameter_has_nowhere_to_apply(self): + with pytest.raises(ValueError, match="nowhere to apply"): + decision_space(fixed={"rcd": 3.0}, bounds={"rcd": (2.0, 5.0)}) diff --git a/tests/unit/calibration/test_runner.py b/tests/unit/calibration/test_runner.py new file mode 100644 index 00000000..b5f88f2a --- /dev/null +++ b/tests/unit/calibration/test_runner.py @@ -0,0 +1,1271 @@ +"""The calibration runner: a small deterministic differential evolution and its artifacts.""" + +import csv +import inspect +import json +import logging +import multiprocessing +import os +import shutil +import subprocess +import sys +from concurrent.futures import ProcessPoolExecutor +from concurrent.futures.process import BrokenProcessPool +from pathlib import Path + +import numpy as np +import pytest + +from rubem import _deps +from rubem.api import Model +from rubem.calibration import runner +from rubem.calibration.objective import Series +from rubem.calibration.parameters import ( + CALIBRATION_PARAMETERS, + FREE_PARAMETERS, + bounds, + parameters_to_vector, +) +from rubem.calibration.runner import ( + EVALUATION_COLUMNS, + OBSERVED_COLUMNS, + STATION_COLUMNS, + CalibrationError, + CalibrationSettings, + _check_observed_series, + _population_size, + _Progress, + _station_ids, + calibrate, +) +from rubem.configuration.model_configuration import ModelConfiguration +from rubem.configuration.model_configuration_file import ModelConfigurationFile +from rubem.configuration.model_configuration_file_v1 import ( + Aggregation, + ModelConfigurationFileV1, +) +from tests.helpers.config import REPO_ROOT +from tests.helpers.synthetic import write_synthetic_dataset + +# SciPy is the optional ``rubem[calibration]`` extra, and an environment +# without it (the frozen one of the byte-exact job, for instance) still +# collects this file; the whole module is skipped there instead of failing to +# import. +scipy_optimize = pytest.importorskip("scipy.optimize") + +TIMESTEPS = 3 + +# Five admissible members, the second one far from the configuration under +# calibration. SciPy replaces the first member with x0, so the configuration +# itself is always part of the initial population. +INIT = np.array( + [ + [1.0, 0.20, 0.20, 0.20, 2.0, 0.20, 0.20, 0.20], + [9.5, 0.90, 0.10, 0.80, 9.0, 0.90, 0.90, 0.90], + [2.0, 0.30, 0.40, 0.40, 3.0, 0.30, 0.30, 0.30], + [6.0, 0.70, 0.25, 0.25, 7.0, 0.70, 0.70, 0.70], + [3.0, 0.40, 0.50, 0.30, 4.0, 0.40, 0.40, 0.40], + ] +) + + +class Dataset: + """The synthetic dataset, its configuration file and the series a plain run wrote.""" + + def __init__(self, tmp_path): + self.config = write_synthetic_dataset(str(tmp_path), timesteps=TIMESTEPS) + self.config_file = tmp_path / "config.json" + self.config_file.write_text(json.dumps(self.config), encoding="utf8") + written = Model.from_config(self.config).run().time_series["arn"][0] + self.observed = tmp_path / "observed.csv" + shutil.copyfile(written, self.observed) + self.run_dir = tmp_path / "calibration" + self.temp_dir = tmp_path / "temp" + self.parameters = ModelConfiguration( + self.config, validate_input=False + ).calibration_parameters.model_dump() + + def settings(self, **overrides): + """The deterministic small search of the tests.""" + defaults = { + "variable": "arn", + "maxiter": 1, + "popsize": 5, + "seed": 1, + "workers": 2, + "temp_dir": str(self.temp_dir), + "init": INIT, + "polish": False, + } + return CalibrationSettings(**{**defaults, **overrides}) + + def calibrate(self, observed=None, **overrides): + return calibrate( + self.config_file, + observed or self.observed, + self.run_dir, + self.settings(**overrides), + ) + + +def read_evaluations(path): + with path.open(encoding="utf-8", newline="") as table: + return list(csv.DictReader(table)) + + +def read_station_table(path): + """Read one of the ``;``-separated per-station tables of a run directory.""" + with path.open(encoding="utf-8", newline="") as table: + return list(csv.DictReader(table, delimiter=";")) + + +@pytest.fixture +def progress_log(caplog): + """Capture the ``rubem.progress`` lines whatever the logging configuration is. + + The command line detaches that logger from the root one so that its lines + are printed verbatim, and a test that ran the command line before this one + may have left it detached; attaching the capture handler to the logger + itself is what makes the capture independent of that. + """ + captured = logging.getLogger("rubem.progress") + handlers, level, propagate = captured.handlers[:], captured.level, captured.propagate + captured.handlers = [caplog.handler] + captured.setLevel(logging.INFO) + captured.propagate = False + caplog.set_level(logging.INFO) + try: + yield caplog + finally: + captured.handlers[:] = handlers + captured.setLevel(level) + captured.propagate = propagate + + +@pytest.fixture +def dataset(tmp_path): + return Dataset(tmp_path) + + +@pytest.fixture(scope="class") +def calibrated(tmp_path_factory): + """One calibration, shared by the assertions that inspect its artifacts.""" + data = Dataset(tmp_path_factory.mktemp("calibration")) + return data, data.calibrate() + + +class TestCalibration: + @pytest.mark.unit + def test_the_configuration_that_produced_the_observations_is_the_optimum(self, calibrated): + data, result = calibrated + + # The observed series is what the configuration itself wrote, so the + # global optimum of the search is the configuration, with NSE 1. + assert result.best_objective == pytest.approx(0.0, abs=1e-9) + assert result.best_nse == pytest.approx(1.0, abs=1e-12) + assert result.best_parameters == pytest.approx(data.parameters, abs=1e-12) + assert set(result.best_parameters) == set(CALIBRATION_PARAMETERS) + + @pytest.mark.unit + def test_every_evaluation_is_a_row_of_the_table(self, calibrated): + data, result = calibrated + rows = read_evaluations(result.evaluations_csv) + + assert list(rows[0]) == list(EVALUATION_COLUMNS) + assert len(rows) == result.evaluations + assert result.best_objective == pytest.approx(min(float(row["objective"]) for row in rows)) + assert {row["id"] for row in rows} == {row["id"] for row in rows if row["id"]} + assert len(rows) >= len(INIT) + # The linear constraint keeps the candidates whose weights do not add up + # away from the workers, so no evaluation is ever rejected as + # inadmissible and none of them failed. + assert [row["error"] for row in rows if row["error"]] == [] + + @pytest.mark.unit + def test_the_evaluations_ran_in_other_processes_and_none_lingers(self, calibrated): + _, result = calibrated + rows = read_evaluations(result.evaluations_csv) + + pids = [int(row["pid"]) for row in rows] + + assert pids, "at least one evaluation was recorded" + assert os.getpid() not in pids + # One task per child: PCRaster keeps its clone and its rasters for the + # life of the process, so no two evaluations may share one. + assert len(set(pids)) == len(pids) + assert multiprocessing.active_children() == [] + + @pytest.mark.unit + def test_the_summary_carries_the_settings_and_the_population_size(self, calibrated): + data, result = calibrated + summary = json.loads(result.result_json.read_text(encoding="utf-8")) + + assert summary["nfev"] == result.evaluations + assert summary["nit"] == result.generations + assert summary["population_size"] == len(INIT) + assert summary["best_objective"] == pytest.approx(result.best_objective) + assert summary["best_nse"] == pytest.approx(result.best_nse) + assert summary["best_parameters"] == pytest.approx(result.best_parameters) + assert summary["settings"]["seed"] == 1 + assert summary["settings"]["workers"] == 2 + assert summary["settings"]["variable"] == "arn" + assert summary["settings"]["maxiter"] == 1 + assert summary["settings"]["polish"] is False + assert summary["settings"]["temp_dir"] == str(data.temp_dir) + assert summary["settings"]["init"] == INIT.tolist() + assert summary["settings"]["allow_blocking_problems"] is False + + @pytest.mark.unit + def test_the_calibrated_configuration_loads_with_the_best_parameters(self, calibrated): + data, result = calibrated + + assert result.calibrated_config.name == "config-calibrated.json" + configuration = ModelConfiguration(result.calibrated_config, validate_input=False) + assert configuration.calibration_parameters.model_dump() == pytest.approx( + result.best_parameters, abs=1e-12 + ) + document = json.loads(result.calibrated_config.read_text(encoding="utf-8")) + assert "CALIBRATION" in document, "a legacy configuration stays legacy" + + @pytest.mark.unit + def test_the_best_row_of_the_table_is_what_was_written_back(self, calibrated): + _, result = calibrated + rows = read_evaluations(result.evaluations_csv) + best = min(rows, key=lambda row: float(row["objective"])) + parameters = {name: float(best[name]) for name in CALIBRATION_PARAMETERS} + + # The optimum SciPy reports and the best record of the run are one and + # the same evaluation, down to the last bit, and it is that evaluation + # the calibrated configuration carries. + assert float(best["objective"]) == result.best_objective + assert float(best["nse"]) == result.best_nse + assert parameters == result.best_parameters + configuration = ModelConfiguration(result.calibrated_config, validate_input=False) + assert configuration.calibration_parameters.model_dump() == pytest.approx( + parameters, abs=1e-12 + ) + + @pytest.mark.unit + def test_the_temporary_directory_is_left_empty(self, calibrated): + data, _ = calibrated + + assert list(data.temp_dir.iterdir()) == [] + + @pytest.mark.unit + def test_the_run_directory_holds_every_artifact(self, calibrated): + _, result = calibrated + + assert result.run_dir.is_dir() + assert result.evaluations_csv.is_file() + assert result.result_json.is_file() + assert result.calibrated_config.is_file() + assert result.stations_csv.is_file() + assert result.best_series.is_file() + assert (result.run_dir / "observed.csv").is_file() + assert result.best_series.name == "best_arn.csv" + assert result.message + assert result.generations == 1 + + @pytest.mark.unit + def test_the_summary_names_every_artifact_it_wrote(self, calibrated): + _, result = calibrated + artifacts = json.loads(result.result_json.read_text(encoding="utf-8"))["artifacts"] + + assert Path(artifacts["evaluations_csv"]) == result.evaluations_csv + assert Path(artifacts["observed_csv"]) == result.run_dir / "observed.csv" + assert Path(artifacts["stations_csv"]) == result.stations_csv + assert Path(artifacts["best_series"]) == result.best_series + assert Path(artifacts["calibrated_config"]) == result.calibrated_config + + @pytest.mark.unit + def test_the_observations_are_summarized_before_the_search(self, calibrated): + _, result = calibrated + rows = read_station_table(result.run_dir / "observed.csv") + + assert list(rows[0]) == list(OBSERVED_COLUMNS) + assert [row["station"] for row in rows] == ["1", "2"] + for row in rows: + # The observed file is the table the configuration itself wrote: + # every compared step is observed and nothing is a gap. + assert row["in_selection"] == "true" + assert row["pairs_in_window"] == str(TIMESTEPS) + assert row["dropped"] == "0" + assert float(row["min"]) <= float(row["mean"]) <= float(row["max"]) + assert float(row["std"]) >= 0.0 + + @pytest.mark.unit + def test_the_stations_of_the_best_candidate_are_a_table_of_their_own(self, calibrated): + _, result = calibrated + rows = read_station_table(result.stations_csv) + + assert list(rows[0]) == list(STATION_COLUMNS) + assert [row["station"] for row in rows] == ["1", "2"] + for row in rows: + # The optimum reproduces the observations exactly. + assert row["in_selection"] == "true" + assert int(row["pairs"]) == TIMESTEPS + assert float(row["nse"]) == pytest.approx(1.0, abs=1e-12) + assert float(row["rmse"]) == pytest.approx(0.0, abs=1e-9) + assert float(row["r"]) == pytest.approx(1.0, abs=1e-9) + assert float(row["mean_observed"]) == pytest.approx(float(row["mean_simulated"])) + assert float(row["std_observed"]) == pytest.approx(float(row["std_simulated"])) + + @pytest.mark.unit + def test_the_series_of_the_best_candidate_is_kept_beside_the_observed_one(self, calibrated): + _, result = calibrated + rows = read_station_table(result.best_series) + + assert list(rows[0]) == [ + "step", + "observed_1", + "simulated_1", + "observed_2", + "simulated_2", + ] + assert [int(row["step"]) for row in rows] == list(range(1, TIMESTEPS + 1)) + for row in rows: + for station in ("1", "2"): + # The best candidate is the configuration that wrote the + # observations, so the two columns of a station are one series. + assert float(row[f"simulated_{station}"]) == pytest.approx( + float(row[f"observed_{station}"]), abs=1e-9 + ) + + @pytest.mark.unit + def test_the_table_of_the_best_candidate_is_read_back_by_the_reader(self, calibrated): + from rubem.calibration.objective import read_series + + _, result = calibrated + table = read_series(result.best_series) + + # It is written in the layout the readers accept, header included. + assert table.steps.tolist() == list(range(1, TIMESTEPS + 1)) + assert sorted(table.stations) == ["observed_1", "observed_2", "simulated_1", "simulated_2"] + + @pytest.mark.unit + def test_every_row_of_the_table_carries_the_moment_it_started(self, calibrated): + import datetime + + _, result = calibrated + rows = read_evaluations(result.evaluations_csv) + + assert "started_at" in EVALUATION_COLUMNS + moments = [datetime.datetime.fromisoformat(row["started_at"]) for row in rows] + assert all(moment.tzinfo is not None for moment in moments) + # The rows are ordered by the candidate they evaluated; the moments are + # what puts them back in the order the evaluations were made in. + assert sorted(moments) != [] and len(moments) == len(rows) + + +class _RecordingExecutor(ProcessPoolExecutor): + """A process pool that remembers how it was built and when it was closed.""" + + built: list["_RecordingExecutor"] = [] + + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.built_with = kwargs + self.shutdowns = 0 + _RecordingExecutor.built.append(self) + + def shutdown(self, *args, **kwargs): + self.shutdowns += 1 + super().shutdown(*args, **kwargs) + + +@pytest.fixture +def broken_search(monkeypatch): + """Record what the search is given, and let it fail without running a model.""" + call = {} + + def explode(function, search_bounds, **kwargs): + call["function"] = function + call["bounds"] = search_bounds + call.update(kwargs) + raise RuntimeError("the search broke") + + _RecordingExecutor.built.clear() + monkeypatch.setattr(scipy_optimize, "differential_evolution", explode) + monkeypatch.setattr(runner, "ProcessPoolExecutor", _RecordingExecutor) + return call + + +class TestProcessTopology: + @pytest.mark.unit + def test_the_pool_is_a_spawned_one_and_is_closed_when_the_search_fails( + self, dataset, broken_search + ): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate() + + (executor,) = _RecordingExecutor.built + assert executor.built_with["max_workers"] == 2 + assert executor.built_with["max_tasks_per_child"] == 1 + assert executor.built_with["mp_context"].get_start_method() == "spawn" + # The pool is closed even though the search raised, and the map of that + # very pool is what the search was told to evaluate with. + assert executor.shutdowns == 1 + assert broken_search["workers"].__self__ is executor + assert multiprocessing.active_children() == [] + + @pytest.mark.unit + def test_the_search_is_given_the_documented_arguments(self, dataset, broken_search): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate() + + assert broken_search["bounds"] == bounds() + assert broken_search["strategy"] == "best1exp" + assert broken_search["maxiter"] == 1 + assert broken_search["popsize"] == 5 + assert broken_search["mutation"] == (0.5, 1.0) + assert broken_search["recombination"] == 0.7 + assert broken_search["updating"] == "deferred" + assert broken_search["polish"] is False + assert broken_search["rng"] == 1 + assert broken_search["init"] is INIT + # The configuration under calibration is the first member of the + # population: SciPy puts x0 in place of the first row of init. + assert broken_search["x0"] == pytest.approx(parameters_to_vector(dataset.parameters)) + + @pytest.mark.unit + def test_the_search_is_constrained_to_weights_that_add_up_to_at_most_one( + self, dataset, broken_search + ): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate() + + constraint = broken_search["constraints"] + coefficients = np.atleast_2d(constraint.A) + expected = [ + [1.0 if name in ("w_1", "w_2") else 0.0 for name in FREE_PARAMETERS], + ] + + assert coefficients.tolist() == expected + assert np.atleast_1d(constraint.ub).tolist() == [1.0] + assert np.all(np.isneginf(np.atleast_1d(constraint.lb))) + + +class TestProgress: + @pytest.mark.unit + def test_the_callback_takes_the_argument_scipy_passes_by_keyword(self): + # SciPy reads the signature of the callback: only a callback whose one + # parameter is called ``intermediate_result`` is given the result of the + # generation, any other is called with the old ``(x, convergence)`` pair. + assert set(inspect.signature(_Progress(Path())).parameters) == {"intermediate_result"} + + @pytest.mark.unit + def test_the_callback_logs_the_generation_and_lets_the_search_go_on( + self, tmp_path, progress_log, capfd + ): + evaluations = tmp_path / "evaluations" + evaluations.mkdir() + (evaluations / "one.json").write_text("{}", encoding="utf8") + callback = _Progress(evaluations) + + stop = callback(scipy_optimize.OptimizeResult(fun=12.5, nit=3)) + + assert stop is False, "the callback never halts the search" + assert callback.generations == 3 + assert "Generation 3" in progress_log.text + assert "12.5" in progress_log.text + assert "1 evaluation" in progress_log.text + # A record without an objective is not the best one; a callback that + # read it as a candidate would end the search it is only watching. + assert "best NSE n/a" in progress_log.text + assert capfd.readouterr().out == "", "the library logs, it does not print" + + @pytest.mark.unit + def test_the_generation_line_carries_the_best_efficiency_of_the_records( + self, tmp_path, progress_log + ): + evaluations = tmp_path / "evaluations" + evaluations.mkdir() + for name, objective_value, nse in (("a", 4.0, 0.25), ("b", 1.0, 0.75)): + (evaluations / f"{name}.json").write_text( + json.dumps({"id": name, "objective": objective_value, "nse": nse, "error": None}), + encoding="utf8", + ) + + _Progress(evaluations)(scipy_optimize.OptimizeResult(fun=1.0, nit=2)) + + # SciPy reports the objective; the efficiency behind it is read from + # the records, which is where the workers write it. + assert "Generation 2" in progress_log.text + assert "best NSE 0.750000" in progress_log.text + assert "2 evaluation(s) recorded" in progress_log.text + + +class TestReproducibility: + @pytest.mark.unit + def test_two_runs_of_one_seed_write_the_same_table(self, dataset, capfd): + first = calibrate( + dataset.config_file, + dataset.observed, + dataset.run_dir / "first", + dataset.settings(), + ) + second = calibrate( + dataset.config_file, + dataset.observed, + dataset.run_dir / "second", + dataset.settings(), + ) + + def candidates(result): + return [ + tuple(row[name] for name in CALIBRATION_PARAMETERS) + for row in read_evaluations(result.evaluations_csv) + ] + + # Same seed, same initial population, same candidates, in the same + # order: the synthetic dataset carries a fixed LDD, so the simulation + # itself is reproducible and the two searches follow the same path. Only + # the identifiers, the process ids and the elapsed times differ. + assert candidates(first) == candidates(second) + assert first.best_parameters == second.best_parameters + assert first.best_objective == second.best_objective + assert first.evaluations == second.evaluations + # The library logs, the command line prints, and the workers are silent: + # the newline the PCRaster framework prints when a worker exits goes to + # the null device, and the warnings of an unvalidated load, which the + # parent already reported, meet a null handler instead of the last + # resort handler of logging. Captured at file-descriptor level, so the + # children are included. + out, err = capfd.readouterr() + assert out == "" + assert "WARNING" not in err and "No handlers could be found" not in err + + +class TestFailures: + @pytest.mark.unit + def test_a_run_whose_evaluations_all_fail_is_reported_with_its_first_error( + self, dataset, tmp_path + ): + # The stations and the steps of this file are the ones the run samples, + # so the parent lets it through; every value is a gap, so every worker + # fails inside the evaluation of the series. + gaps = tmp_path / "gaps.csv" + gaps.write_text( + "0;1;2\n" + "".join(f"{step};-9999;-9999\n" for step in range(1, TIMESTEPS + 1)), + encoding="utf8", + ) + + with pytest.raises(CalibrationError, match="Nash-Sutcliffe") as failure: + dataset.calibrate(observed=gaps, maxiter=1) + + # The search ran: what is reported is the summary of the finished run, + # carrying the first error a worker recorded, not a refusal of the parent. + assert "Every evaluation of the calibration failed" in str(failure.value) + + rows = read_evaluations(dataset.run_dir / "evaluations.csv") + assert rows, "the table of the evaluations is written even when every one failed" + assert all(row["error"] for row in rows) + assert all(row["nse"] == "" for row in rows) + + @pytest.mark.unit + def test_a_variable_the_model_does_not_write_is_refused(self, dataset): + with pytest.raises(CalibrationError, match="not an output variable"): + dataset.calibrate(variable="flow") + + @pytest.mark.unit + def test_a_run_directory_of_an_earlier_calibration_is_refused(self, dataset): + evaluations = dataset.run_dir / "evaluations" + evaluations.mkdir(parents=True) + (evaluations / "deadbeef.json").write_text("{}", encoding="utf8") + + with pytest.raises(CalibrationError, match="already holds 1 evaluation"): + dataset.calibrate() + + @pytest.mark.unit + def test_a_negative_spin_up_is_refused(self, dataset): + with pytest.raises(CalibrationError, match="must not be negative"): + dataset.calibrate(spinup_steps=-1) + + @pytest.mark.unit + def test_a_search_without_a_worker_is_refused(self, dataset): + with pytest.raises(CalibrationError, match="at least one worker"): + dataset.calibrate(workers=0) + + @pytest.mark.unit + def test_without_scipy_the_optional_extra_is_named(self, dataset, monkeypatch): + monkeypatch.setattr(_deps, "missing_calibration_deps", lambda: ["scipy"]) + + with pytest.raises(CalibrationError, match=r"rubem\[calibration\]"): + dataset.calibrate() + + assert not dataset.run_dir.exists(), "nothing is written when the search cannot start" + + +@pytest.fixture +def no_pool(monkeypatch): + """Let a test prove that no worker process was ever started.""" + _RecordingExecutor.built.clear() + monkeypatch.setattr(runner, "ProcessPoolExecutor", _RecordingExecutor) + return _RecordingExecutor.built + + +class TestObservedSeries: + @pytest.mark.unit + def test_an_observed_series_of_foreign_stations_is_refused_before_the_search( + self, dataset, tmp_path, no_pool + ): + foreign = tmp_path / "foreign.csv" + foreign.write_text("0;A;B\n1;1.0;2.0\n2;2.0;3.0\n3;3.0;4.0\n", encoding="utf8") + + with pytest.raises(CalibrationError, match="no station in common") as failure: + dataset.calibrate(observed=foreign, maxiter=1) + + # Both sets are named, and the refusal comes before anything is written + # or any worker is started: a wrong observed file costs no model run. + assert "A, B" in str(failure.value) + assert "1, 2" in str(failure.value) + assert not dataset.run_dir.exists() + assert no_pool == [] + assert multiprocessing.active_children() == [] + + @pytest.mark.unit + @pytest.mark.parametrize( + ("steps", "overrides"), + [ + ([TIMESTEPS + 1, TIMESTEPS + 2], {}), + ([-1, 0], {}), + ([1, TIMESTEPS], {"spinup_steps": TIMESTEPS}), + ], + ids=["after the window", "before the window", "inside the spin-up"], + ) + def test_an_observed_series_outside_the_compared_steps_is_refused( + self, dataset, tmp_path, no_pool, steps, overrides + ): + elsewhere = tmp_path / "elsewhere.csv" + elsewhere.write_text( + "0;1;2\n" + "".join(f"{step};1.0;2.0\n" for step in steps), encoding="utf8" + ) + + with pytest.raises(CalibrationError, match="no time step in common") as failure: + dataset.calibrate(observed=elsewhere, **overrides) + + message = str(failure.value) + assert f"{min(steps)} to {max(steps)}" in message + assert f"1 to {TIMESTEPS}" in message + assert f"{overrides.get('spinup_steps', 0)} spin-up" in message + assert not dataset.run_dir.exists() + assert no_pool == [] + + @pytest.mark.unit + def test_an_observed_station_the_run_does_not_sample_is_warned_about_and_ignored( + self, dataset, caplog + ): + rows = dataset.observed.read_text(encoding="utf8").splitlines() + extra = dataset.observed.with_name("extra.csv") + extra.write_text( + "".join(f"{row};{'3' if number == 0 else '1.0'}\n" for number, row in enumerate(rows)), + encoding="utf8", + ) + + exact = calibrate( + dataset.config_file, dataset.observed, dataset.run_dir / "exact", dataset.settings() + ) + with caplog.at_level(logging.WARNING, logger="rubem.calibration.runner"): + with_extra = calibrate( + dataset.config_file, extra, dataset.run_dir / "extra", dataset.settings() + ) + + # The foreign column is named once and then ignored by the alignment, so + # the calibration is the one the exact file produces. + assert "does not sample (3)" in caplog.text + assert with_extra.best_parameters == exact.best_parameters + assert with_extra.best_objective == exact.best_objective + assert with_extra.best_nse == exact.best_nse + assert with_extra.evaluations == exact.evaluations + + @pytest.mark.unit + def test_a_negative_observation_is_counted_as_a_gap_before_the_search( + self, dataset, broken_search, progress_log + ): + rows = dataset.observed.read_text(encoding="utf8").splitlines() + first = rows[1].split(";") + # The gauge of station 1 has no reading at the first step, written as + # the negative marker an observation table marks a gap with. + first[1] = "-1.5" + gapped = dataset.observed.with_name("gapped.csv") + gapped.write_text("\n".join([rows[0], ";".join(first), *rows[2:]]) + "\n", encoding="utf8") + + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate(observed=gapped) + + summarized = { + row["station"]: row for row in read_station_table(dataset.run_dir / "observed.csv") + } + # The negative value never reaches the statistics, and both the table + # and the line the run prints count it as a gap of that station alone. + assert summarized["1"]["dropped"] == "1" + assert summarized["1"]["pairs_in_window"] == str(TIMESTEPS - 1) + assert float(summarized["1"]["min"]) > 0.0 + assert summarized["2"]["dropped"] == "0" + assert summarized["2"]["pairs_in_window"] == str(TIMESTEPS) + assert ( + f"Station 1: {TIMESTEPS - 1} observed value(s) on the compared steps, " + "1 dropped as gaps." in progress_log.text + ) + + @pytest.mark.unit + def test_the_stations_are_the_values_of_the_sample_locations_raster(self, dataset): + samples = ModelConfiguration( + dataset.config, validate_input=False + ).raster_files.sample_locations + + # The no-data cells of the nominal map are not stations, and the ids are + # the strings the model writes in the header of its table. + assert _station_ids(samples) == {"1", "2"} + + @pytest.mark.unit + def test_a_background_of_zeroes_is_not_a_station(self, tmp_path): + from osgeo import gdal + + gdal.UseExceptions() + gdal.AllRegister() + raster = tmp_path / "samples-with-zeroes.tif" + handle = gdal.GetDriverByName("GTiff").Create(str(raster), 3, 3, 1, gdal.GDT_Int32) + try: + values = np.zeros((3, 3), dtype=np.int32) + values[0, 0], values[2, 2] = 1, 2 + handle.GetRasterBand(1).WriteArray(values) + handle.FlushCache() + finally: + handle = None + + # A raster may mark its background with zeroes and declare no no-data + # value at all; the validation of the inputs counts neither the zeroes + # nor the missing cells as a sample, and neither does this. + assert _station_ids(raster) == {"1", "2"} + + @pytest.mark.unit + def test_the_stations_of_a_zones_aggregation_are_not_checked(self): + # The model remaps the zone ids to 1..N at run time, so the parent has + # no station id to compare and the check only covers the steps. + observed = Series(steps=np.array([1, 2, 3]), stations={"A": np.zeros(3)}) + + assert ( + _check_observed_series( + observed, + first_step=1, + last_step=3, + spinup_steps=0, + aggregation=Aggregation.ZONES, + station_ids={"1", "2"}, + ) + is None + ) + + +class TestSettings: + @pytest.mark.unit + def test_the_default_decision_vector_is_the_configuration_under_calibration(self, dataset): + vector = parameters_to_vector(dataset.parameters) + + assert vector == pytest.approx([4.5, 0.5, 0.333, 0.333, 5.0, 0.5, 0.5, 0.5]) + + @pytest.mark.unit + def test_the_defaults_are_the_documented_ones(self): + settings = CalibrationSettings() + + assert (settings.variable, settings.spinup_steps) == ("arn", 0) + assert (settings.maxiter, settings.popsize) == (100, 15) + assert settings.init == "sobol" + assert settings.strategy == "best1exp" + assert settings.mutation == (0.5, 1.0) + assert settings.recombination == 0.7 + assert settings.polish is False + assert settings.workers is None + assert settings.allow_blocking_problems is False + + @pytest.mark.unit + @pytest.mark.parametrize("allowed", [False, True]) + def test_allow_blocking_problems_reaches_the_configuration_loader( + self, dataset, monkeypatch, allowed + ): + """The single load of a calibration carries the option to ``Model.from_file``. + + The workers never revalidate, so this call is the only place the option + can have an effect, and the search itself is irrelevant to the question. + """ + seen = {} + original = Model.from_file + + def recording_from_file(path, **keywords): + seen.update(keywords) + return original(path, **keywords) + + def stop(*_args, **_keywords): + raise RuntimeError("the load is all this test exercises") + + monkeypatch.setattr(Model, "from_file", staticmethod(recording_from_file)) + # The search that follows the load is irrelevant here, so the next step + # of the calibration ends it. + monkeypatch.setattr(runner, "read_series", stop) + + with pytest.raises(RuntimeError, match="the load is all this test exercises"): + runner.calibrate( + dataset.config_file, + dataset.observed, + dataset.run_dir / "allowed", + CalibrationSettings(allow_blocking_problems=allowed), + ) + + assert seen["allow_blocking_problems"] is allowed + assert seen["validate_input"] is True + + +@pytest.fixture(scope="class") +def calibrated_v1(tmp_path_factory): + """One calibration of the same dataset written as a format 1.0 configuration.""" + tmp_path = tmp_path_factory.mktemp("calibration-v1") + data = Dataset(tmp_path) + legacy = ModelConfigurationFile.model_validate(data.config) + data.config_file = tmp_path / "basin.json" + data.config_file.write_text( + json.dumps(ModelConfigurationFileV1.from_legacy(legacy).to_dict()), encoding="utf8" + ) + return data, data.calibrate() + + +class TestFormatV1Configuration: + @pytest.mark.unit + def test_a_format_1_0_configuration_is_calibrated_and_written_back_in_its_format( + self, calibrated_v1 + ): + data, result = calibrated_v1 + + assert result.best_objective == pytest.approx(0.0, abs=1e-9) + assert result.best_parameters == pytest.approx(data.parameters, abs=1e-12) + assert result.calibrated_config.name == "basin-calibrated.json" + document = json.loads(result.calibrated_config.read_text(encoding="utf-8")) + assert document["version"] == "1.0" + assert document["model_calibration_parameters"]["b"] == pytest.approx( + result.best_parameters["beta"], abs=1e-12 + ) + configuration = ModelConfiguration(result.calibrated_config, validate_input=False) + assert configuration.calibration_parameters.model_dump() == pytest.approx( + result.best_parameters, abs=1e-12 + ) + + +class TestEvaluationBudget: + @pytest.mark.unit + def test_a_sobol_population_is_rounded_up_to_a_power_of_two(self): + # The documented arithmetic: max(5, popsize * 8) members, rounded up to + # the next power of two, because a Sobol' sequence is balanced only over + # a power-of-two sample. + assert _population_size(CalibrationSettings()) == 128 + assert _population_size(CalibrationSettings(popsize=1)) == 8 + assert _population_size(CalibrationSettings(popsize=16)) == 128 + + @pytest.mark.unit + def test_another_initialization_uses_the_members_it_asks_for(self): + assert _population_size(CalibrationSettings(init="latinhypercube", popsize=3)) == 24 + assert _population_size(CalibrationSettings(init="random", popsize=15)) == 120 + assert _population_size(CalibrationSettings(init=INIT)) == len(INIT) + + @pytest.mark.unit + def test_a_smaller_decision_space_needs_fewer_members(self): + # One fixed parameter takes a dimension out of the search, and the + # number of members follows it. + assert _population_size(CalibrationSettings(popsize=15), 7) == 128 + assert _population_size(CalibrationSettings(popsize=15, init="latinhypercube"), 7) == 105 + assert _population_size(CalibrationSettings(popsize=1, init="random"), 1) == 5 + + @pytest.mark.unit + @pytest.mark.parametrize("popsize", [1, 5, 15, 16]) + @pytest.mark.parametrize("init", ["sobol", "halton", "latinhypercube", "random"]) + @pytest.mark.parametrize("dimension", [7, 8]) + def test_the_announced_population_is_the_one_scipy_builds(self, popsize, init, dimension): + # The announced number is what the command line reports as the budget, + # so it has to be SciPy's own arithmetic and not an approximation of it: + # only ``sobol`` rounds the population up to a power of two, and the + # number of free parameters is a factor of it. + optimum = scipy_optimize.differential_evolution( + lambda vector: float(np.sum(vector)), + bounds()[:dimension], + init=init, + popsize=popsize, + maxiter=0, + polish=False, + rng=0, + ) + + assert len(optimum.population) == _population_size( + CalibrationSettings(popsize=popsize, init=init), dimension + ) + + @pytest.mark.unit + def test_the_summary_reports_the_population_the_search_really_built(self, dataset, monkeypatch): + # A named initialization resizes the population, so the summary must + # report the population of the result and not the number this module + # announced for the budget. + def stub(function, search_bounds, **kwargs): + population = np.asarray(kwargs["init"], dtype=np.float64) + energies = np.asarray([function(member) for member in population[:2]]) + enlarged = np.repeat(population, 2, axis=0) + return scipy_optimize.OptimizeResult( + x=population[0], + fun=float(energies[0]), + nfev=len(energies), + nit=1, + success=True, + message="stub", + population=enlarged, + population_energies=np.full(len(enlarged), np.inf), + ) + + monkeypatch.setattr(scipy_optimize, "differential_evolution", stub) + + result = dataset.calibrate() + summary = json.loads(result.result_json.read_text(encoding="utf-8")) + + assert summary["population_size"] == 2 * len(INIT) + assert _population_size(dataset.settings()) == len(INIT) + + +class TestOptionalDependency: + @pytest.mark.unit + def test_importing_the_package_does_not_import_scipy(self, tmp_path): + # An installation without the rubem[calibration] extra must still be + # able to import the package, so SciPy is imported inside the functions + # that need it. A fresh interpreter is the only honest probe: the test + # session has SciPy loaded already. + probe = tmp_path / "probe.py" + probe.write_text( + "import sys\n" + "import rubem.calibration.runner\n" + "import rubem.calibration._worker\n" + "import rubem.calibration.parameters\n" + "print('scipy' in sys.modules)\n", + encoding="utf8", + ) + + environment = {**os.environ, "PYTHONPATH": REPO_ROOT} + completed = subprocess.run( + [sys.executable, str(probe)], + capture_output=True, + text=True, + check=True, + cwd=REPO_ROOT, + env=environment, + ) + + assert completed.stdout.strip() == "False" + + +@pytest.fixture(scope="class") +def calibrated_with_a_fixed_parameter(tmp_path_factory): + """One calibration in which ``x`` is pinned and only seven parameters are searched.""" + data = Dataset(tmp_path_factory.mktemp("calibration-fixed")) + free = np.delete(INIT, FREE_PARAMETERS.index("x"), axis=1) + return data, data.calibrate(fixed={"x": 0.0}, init=free, stations=("1", "2")) + + +class TestFixedParameters: + @pytest.mark.unit + def test_the_fixed_parameter_is_not_searched_and_is_in_every_candidate( + self, calibrated_with_a_fixed_parameter + ): + _, result = calibrated_with_a_fixed_parameter + rows = read_evaluations(result.evaluations_csv) + + assert set(result.best_parameters) == set(CALIBRATION_PARAMETERS) + assert result.best_parameters["x"] == 0.0 + assert rows, "the search ran" + assert [row["error"] for row in rows if row["error"]] == [] + # Every candidate the search proposed carries the pinned value, and the + # column is still in the table: the parameter left the decision vector, + # not the model. + assert {row["x"] for row in rows} == {"0.0"} + + @pytest.mark.unit + def test_the_summary_records_the_fixed_value_the_bounds_and_the_stations( + self, calibrated_with_a_fixed_parameter + ): + _, result = calibrated_with_a_fixed_parameter + settings = json.loads(result.result_json.read_text(encoding="utf-8"))["settings"] + + assert settings["fixed"] == {"x": 0.0} + assert settings["stations"] == ["1", "2"] + assert list(settings["bounds"]) == [name for name in FREE_PARAMETERS if name != "x"] + assert settings["bounds"]["rcd"] == [1.0, 10.0] + + @pytest.mark.unit + def test_the_calibrated_configuration_carries_the_fixed_value( + self, calibrated_with_a_fixed_parameter + ): + _, result = calibrated_with_a_fixed_parameter + configuration = ModelConfiguration(result.calibrated_config, validate_input=False) + + assert configuration.calibration_parameters.model_dump()["x"] == 0.0 + + @pytest.mark.unit + def test_a_fixed_weight_replaces_the_constraint_with_a_narrowed_bound( + self, dataset, broken_search + ): + free = np.delete(INIT, FREE_PARAMETERS.index("w_1"), axis=1) + + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate(fixed={"w_1": 0.7}, init=free) + + searched = [name for name in FREE_PARAMETERS if name != "w_1"] + # One free weight is a bound, not a constraint: the argument is left out + # instead of being handed over as ``None``. + assert "constraints" not in broken_search + assert len(broken_search["bounds"]) == len(searched) == 7 + assert broken_search["bounds"][searched.index("w_2")] == (0.0, pytest.approx(0.3)) + assert len(broken_search["x0"]) == 7 + + @pytest.mark.unit + def test_a_decision_space_the_settings_cannot_describe_is_refused(self, dataset, no_pool): + with pytest.raises(ValueError, match="cannot be fixed"): + dataset.calibrate(fixed={"x": 5.0}) + + assert not dataset.run_dir.exists() + assert no_pool == [] + + @pytest.mark.unit + def test_an_initial_population_of_the_wrong_width_is_refused(self, dataset, no_pool): + # ``INIT`` has one column per parameter of the whole search, and the + # fixed one has left the decision vector. + with pytest.raises(ValueError, match="one column per free parameter") as failure: + dataset.calibrate(fixed={"x": 0.0}) + + assert "(5, 8)" in str(failure.value) + assert "7 free parameter(s)" in str(failure.value) + assert not dataset.run_dir.exists() + assert no_pool == [] + + @pytest.mark.unit + def test_an_initial_population_that_is_not_a_table_of_members_is_refused( + self, dataset, no_pool + ): + with pytest.raises(ValueError, match="two-dimensional array"): + dataset.calibrate(init=INIT[0]) + + assert not dataset.run_dir.exists() + assert no_pool == [] + + +class TestNarrowedBounds: + @pytest.mark.unit + def test_the_search_is_given_the_narrowed_bounds(self, dataset, broken_search): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate(bounds={"rcd": (2.0, 5.0)}) + + expected = bounds() + expected[FREE_PARAMETERS.index("rcd")] = (2.0, 5.0) + assert broken_search["bounds"] == expected + + @pytest.mark.unit + def test_a_starting_point_outside_a_narrowed_bound_is_moved_onto_it( + self, dataset, broken_search, caplog + ): + with caplog.at_level(logging.WARNING, logger="rubem.calibration.runner"): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate(bounds={"alpha": (5.0, 8.0)}) + + # The configuration has alpha = 4.5, which this run does not search; + # SciPy refuses a starting point outside the bounds, so it is moved onto + # the bound and the move is reported instead of ending the run. + assert broken_search["x0"][FREE_PARAMETERS.index("alpha")] == 5.0 + assert "alpha" in caplog.text + assert "4.5" in caplog.text + + @pytest.mark.unit + def test_a_bound_that_does_not_narrow_the_settings_range_is_refused(self, dataset, no_pool): + with pytest.raises(ValueError, match="not a narrower range"): + dataset.calibrate(bounds={"rcd": (0.5, 12.0)}) + + assert not dataset.run_dir.exists() + assert no_pool == [] + + +class TestStationSelection: + @pytest.mark.unit + @pytest.mark.parametrize( + "stations", + [(), ("1", ""), ("1", 2), "12"], + ids=["empty", "a blank id", "not a string", "a bare string"], + ) + def test_a_selection_that_names_no_station_is_refused(self, dataset, no_pool, stations): + with pytest.raises(CalibrationError): + dataset.calibrate(stations=stations) + + assert not dataset.run_dir.exists() + assert no_pool == [] + + @pytest.mark.unit + def test_the_selected_stations_are_logged_before_the_search( + self, dataset, broken_search, caplog + ): + with caplog.at_level(logging.INFO, logger="rubem.calibration.runner"): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate(stations=("1",)) + + assert "the station(s) 1" in caplog.text + + +def _record_then_raise(error): + """A search that records one evaluation and then fails, without a worker. + + The records are what an interrupted calibration has to leave behind, and + the failure has to happen where the real one does, inside the search. The + evaluations directory is read out of the objective the search was handed, + which is where the parent put it. + """ + + def stub(function, search_bounds, **kwargs): + evaluations_dir = Path(function.keywords["context"].evaluations_dir) + record = { + "id": "deadbeef", + "pid": 4321, + "started_at": "2000-01-01T00:00:00+00:00", + "parameters": dict.fromkeys(CALIBRATION_PARAMETERS, 0.5), + "nse": 0.25, + "station_nse": {"1": 0.25}, + "station_metrics": {}, + "objective": 5625000.0, + "elapsed_seconds": 1.5, + "error": None, + } + (evaluations_dir / "deadbeef.json").write_text(json.dumps(record), encoding="utf8") + raise error + + return stub + + +class TestInterruptedSearch: + @pytest.mark.unit + def test_an_interrupted_search_still_writes_the_evaluations_it_made(self, dataset, monkeypatch): + monkeypatch.setattr( + scipy_optimize, "differential_evolution", _record_then_raise(KeyboardInterrupt()) + ) + + with pytest.raises(KeyboardInterrupt): + dataset.calibrate() + + # The evaluations are hours of model runs on a real basin; a + # calibration that has to be stopped is read from them like any other. + rows = read_evaluations(dataset.run_dir / "evaluations.csv") + assert [row["id"] for row in rows] == ["deadbeef"] + assert rows[0]["started_at"] == "2000-01-01T00:00:00+00:00" + assert rows[0]["nse"] == repr(0.25) + + @pytest.mark.unit + def test_a_crashed_search_still_writes_the_evaluations_it_made(self, dataset, monkeypatch): + monkeypatch.setattr( + scipy_optimize, "differential_evolution", _record_then_raise(RuntimeError("boom")) + ) + + with pytest.raises(RuntimeError, match="boom"): + dataset.calibrate() + + assert read_evaluations(dataset.run_dir / "evaluations.csv") + + @pytest.mark.unit + def test_a_search_that_recorded_nothing_leaves_no_table(self, dataset, broken_search): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate() + + assert not (dataset.run_dir / "evaluations.csv").exists() + + @pytest.mark.unit + def test_a_worker_killed_by_the_system_is_reported_with_the_advice(self, dataset, monkeypatch): + monkeypatch.setattr( + scipy_optimize, + "differential_evolution", + _record_then_raise(BrokenProcessPool("A process in the process pool was terminated")), + ) + + with pytest.raises(CalibrationError, match="--workers") as failure: + dataset.calibrate() + + message = str(failure.value) + assert "died before finishing" in message + assert "memory" in message + # The cause is kept, so the traceback still says what the pool reported. + assert isinstance(failure.value.__cause__, BrokenProcessPool) + assert read_evaluations(dataset.run_dir / "evaluations.csv") + + +class TestProgressReporting: + @pytest.mark.unit + def test_the_budget_the_coverage_the_generations_and_the_summary_are_reported( + self, dataset, progress_log + ): + result = dataset.calibrate() + + text = progress_log.text + # The budget, before anything is spent on it. + assert "8 free parameter(s)" in text + assert f"{len(INIT)} population member(s)" in text + assert "2 worker process(es)" in text + # The coverage of the compared window, per station. + assert f"compares {TIMESTEPS} time step(s)" in text + assert f"covers {TIMESTEPS} of them" in text + assert f"Station 1: {TIMESTEPS} observed value(s)" in text + assert "0 dropped as gaps" in text + # One line per generation, and the closing summary. + assert "Generation 1:" in text + assert "best NSE" in text + assert "Calibration finished after" in text + assert f"{result.evaluations} evaluation(s)" in text + + @pytest.mark.unit + def test_the_progress_is_logged_and_never_printed(self, dataset, progress_log, capfd): + dataset.calibrate() + + assert "Generation 1:" in progress_log.text + out, _ = capfd.readouterr() + assert out == "", "the library logs, the command line prints" + + +class TestStationsOfTheObjective: + @pytest.mark.unit + def test_a_station_left_out_of_the_objective_is_measured_and_marked(self, dataset): + result = dataset.calibrate(stations=("1",)) + + observed = { + row["station"]: row for row in read_station_table(result.run_dir / "observed.csv") + } + stations = {row["station"]: row for row in read_station_table(result.stations_csv)} + + assert observed["1"]["in_selection"] == "true" + assert observed["2"]["in_selection"] == "false" + # The station outside the objective is measured all the same: that is + # the validation half of a calibration/validation split. + assert stations["1"]["in_selection"] == "true" + assert stations["2"]["in_selection"] == "false" + assert float(stations["2"]["nse"]) == pytest.approx(1.0, abs=1e-12) + assert int(stations["2"]["pairs"]) == TIMESTEPS + assert json.loads(result.result_json.read_text(encoding="utf-8"))["settings"][ + "stations" + ] == ["1"] + + @pytest.mark.unit + def test_a_selection_of_stations_the_run_does_not_have_is_refused_before_the_search( + self, dataset, no_pool + ): + with pytest.raises(CalibrationError, match="none of the station") as failure: + dataset.calibrate(stations=("7", "8")) + + # A misspelt id would otherwise fail every evaluation of the run and + # spend the whole budget finding out. + assert "7, 8" in str(failure.value) + assert "1, 2" in str(failure.value) + assert not dataset.run_dir.exists() + assert no_pool == [] + + @pytest.mark.unit + def test_a_selection_that_names_one_unknown_station_warns_and_keeps_the_rest( + self, dataset, broken_search, caplog + ): + with caplog.at_level(logging.WARNING, logger="rubem.calibration.runner"): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate(stations=("1", "7")) + + assert "does not have (7)" in caplog.text diff --git a/tests/unit/calibration/test_worker.py b/tests/unit/calibration/test_worker.py new file mode 100644 index 00000000..35948e6c --- /dev/null +++ b/tests/unit/calibration/test_worker.py @@ -0,0 +1,383 @@ +"""One evaluation of a candidate: the model run, the record and the failure modes.""" + +import datetime +import json +import os +import shutil +from pathlib import Path + +import pytest + +from rubem.api import Model +from rubem.calibration._worker import ( + EvaluationContext, + _derived_document, + evaluate, + simulate_best, +) +from rubem.calibration.objective import INADMISSIBLE_OBJECTIVE, read_series +from rubem.calibration.parameters import ( + FREE_PARAMETERS, + parameters_to_vector, + vector_to_parameters, +) +from rubem.configuration.model_configuration import ModelConfiguration +from rubem.configuration.model_configuration_file_v1 import ( + VARIABLE_IDS, + ModelConfigurationFileV1, +) +from tests.helpers.synthetic import write_synthetic_dataset + +TIMESTEPS = 3 + + +class Dataset: + """The synthetic dataset, its configuration and the series a plain run wrote.""" + + def __init__(self, tmp_path): + self.config = write_synthetic_dataset(str(tmp_path), timesteps=TIMESTEPS) + self.observed_path = Model.from_config(self.config).run().time_series["arn"][0] + self.configuration = ModelConfiguration(self.config, validate_input=False) + self.document = ModelConfigurationFileV1.from_legacy(self.configuration.file).to_dict() + self.vector = parameters_to_vector(self.configuration.calibration_parameters.model_dump()) + self.temp_dir = tmp_path / "temp" + self.evaluations_dir = tmp_path / "evaluations" + self.temp_dir.mkdir() + self.evaluations_dir.mkdir() + + def context(self, observed_path=None, variable="arn", spinup_steps=0, stations=None): + """An evaluation context reading ``observed_path``, the plain run by default.""" + return EvaluationContext( + document=self.document, + base_dir=self.configuration.base_dir, + variable=variable, + observed=read_series(observed_path or self.observed_path), + spinup_steps=spinup_steps, + temp_dir=str(self.temp_dir), + evaluations_dir=str(self.evaluations_dir), + validate_input=False, + stations=stations, + ) + + def records(self): + """Every record written so far, without the files still being staged.""" + return [ + json.loads(path.read_text(encoding="utf-8")) + for path in sorted(self.evaluations_dir.glob("*.json")) + if not path.name.startswith(".") + ] + + def staged(self): + """The staging files left behind by an interrupted write, if any.""" + return [path.name for path in self.evaluations_dir.iterdir() if ".tmp" in path.name] + + +@pytest.fixture +def dataset(tmp_path): + return Dataset(tmp_path) + + +class TestSuccessfulEvaluation: + @pytest.mark.unit + def test_the_configuration_of_the_observations_scores_a_perfect_objective(self, dataset): + # The observed series is what the configuration itself produced, so its + # own parameters reproduce it exactly: NSE 1, objective 0. + value = evaluate(dataset.vector, dataset.context()) + + assert value == pytest.approx(0.0, abs=1e-9) + (record,) = dataset.records() + assert record["error"] is None + assert record["nse"] == pytest.approx(1.0, abs=1e-12) + assert record["objective"] == pytest.approx(0.0, abs=1e-9) + assert record["station_nse"] == {"1": pytest.approx(1.0), "2": pytest.approx(1.0)} + assert record["parameters"]["w_3"] == pytest.approx(0.334, abs=1e-12) + assert record["pid"] == os.getpid() + assert record["elapsed_seconds"] >= 0.0 + + @pytest.mark.unit + def test_the_temporary_output_directory_is_removed_after_the_run(self, dataset): + evaluate(dataset.vector, dataset.context()) + + assert list(dataset.temp_dir.iterdir()) == [] + + @pytest.mark.unit + def test_the_record_is_installed_by_a_rename_and_leaves_no_staging_file(self, dataset): + evaluate(dataset.vector, dataset.context()) + + assert dataset.staged() == [] + assert len(dataset.records()) == 1 + + @pytest.mark.unit + def test_a_different_candidate_is_evaluated_and_scores_worse(self, dataset): + vector = dataset.vector.copy() + vector[FREE_PARAMETERS.index("alpha_gw")] = 0.01 + + value = evaluate(vector, dataset.context()) + + assert 0.0 < value < INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert record["error"] is None + assert record["nse"] < 1.0 + + +class TestDerivedConfiguration: + @pytest.mark.unit + def test_it_disables_every_raster_series_and_keeps_one_time_series(self, dataset): + document = _derived_document( + dataset.context(), vector_to_parameters(dataset.vector), "/output" + ) + output = document["model_simulation_output"] + rasters = output["raster_series"] + samples = output["time_series_samples"] + + assert output["dir_path"] == "/output" + assert rasters["formats"] == [] + assert [name for name in VARIABLE_IDS if rasters.get(name)] == [] + assert [name for name in VARIABLE_IDS if samples.get(name)] == ["arn"] + assert samples["formats"] == ["CSV"] + assert samples["aggregation"] == "point" + # What the worker hands to the model has to be a configuration again. + ModelConfigurationFileV1.model_validate(document) + + @pytest.mark.unit + def test_the_run_writes_nothing_but_the_table_of_the_calibrated_variable( + self, dataset, monkeypatch + ): + # The output directory is listed on its way out, so that the listing is + # of the directory the evaluation really wrote and the removal still + # happens exactly as it does without the test. + listed = [] + removal = shutil.rmtree + + def spy(path, **kwargs): + listed.extend(sorted(entry.name for entry in Path(path).rglob("*"))) + removal(path, **kwargs) + + monkeypatch.setattr(shutil, "rmtree", spy) + + evaluate(dataset.vector, dataset.context()) + + assert listed == ["metadata.json", "tss_arn.csv"] + + +class TestRejectedCandidate: + @pytest.mark.unit + def test_an_inadmissible_candidate_is_recorded_without_a_model_run(self, dataset): + vector = dataset.vector.copy() + vector[FREE_PARAMETERS.index("w_1")] = 0.6 + vector[FREE_PARAMETERS.index("w_2")] = 0.6 + + value = evaluate(vector, dataset.context()) + + assert value == INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert record["error"] == "inadmissible" + assert record["nse"] is None + assert record["objective"] == INADMISSIBLE_OBJECTIVE + assert list(dataset.temp_dir.iterdir()) == [], "no output directory was created" + + +class TestRecord: + @pytest.mark.unit + def test_the_record_is_staged_in_the_directory_it_is_renamed_into(self, dataset, monkeypatch): + # A rename is atomic only inside one filesystem, so the staging file has + # to be a sibling of the record and not a file in the system temporary + # directory, which may be a different mount. + renames = [] + rename = Path.replace + + def spy(self, target): + renames.append((Path(self), Path(target))) + return rename(self, target) + + monkeypatch.setattr(Path, "replace", spy) + + evaluate(dataset.vector, dataset.context()) + + installed = [pair for pair in renames if pair[1].parent == dataset.evaluations_dir] + assert len(installed) == 1, "the record was installed by exactly one rename" + temporary, target = installed[0] + assert temporary.parent == dataset.evaluations_dir + assert temporary.name.startswith(".") + assert target.name.endswith(".json") + + +class TestFailedEvaluation: + @pytest.mark.unit + def test_observations_of_another_station_are_recorded_as_the_error(self, dataset, tmp_path): + foreign = tmp_path / "foreign.csv" + foreign.write_text("0;A;B\n1;1.0;2.0\n2;2.0;3.0\n3;3.0;4.0\n", encoding="utf8") + + value = evaluate(dataset.vector, dataset.context(observed_path=foreign)) + + assert value == INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert "no station in common" in record["error"] + assert record["error"].startswith("ValueError: ") + assert record["nse"] is None + assert list(dataset.temp_dir.iterdir()) == [], "the output directory was removed anyway" + + @pytest.mark.unit + def test_a_variable_the_run_does_not_write_is_recorded_as_the_error(self, dataset): + value = evaluate(dataset.vector, dataset.context(variable="nope")) + + assert value == INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert record["error"] + assert list(dataset.temp_dir.iterdir()) == [] + + +class TestRecordedMetrics: + @pytest.mark.unit + def test_the_record_carries_the_moment_the_evaluation_started(self, dataset): + before = datetime.datetime.now(datetime.UTC) + + evaluate(dataset.vector, dataset.context()) + + (record,) = dataset.records() + started_at = datetime.datetime.fromisoformat(record["started_at"]) + assert started_at.tzinfo is not None, "the moment carries its time zone" + assert started_at.utcoffset() == datetime.timedelta(0) + # The wall-clock moment is what puts the records back in the order they + # were made in, so it has to be the moment of this evaluation. + assert before <= started_at <= datetime.datetime.now(datetime.UTC) + + @pytest.mark.unit + def test_the_record_carries_the_goodness_of_fit_of_every_station(self, dataset): + evaluate(dataset.vector, dataset.context()) + + (record,) = dataset.records() + metrics = record["station_metrics"] + + assert sorted(metrics) == ["1", "2"] + for measured in metrics.values(): + assert sorted(measured) == [ + "mean_observed", + "mean_simulated", + "nse", + "pairs", + "r", + "rmse", + "std_observed", + "std_simulated", + ] + # The observations are what this very configuration wrote, so the + # simulation reproduces them exactly. + assert measured["pairs"] == TIMESTEPS + assert measured["nse"] == pytest.approx(1.0, abs=1e-12) + assert measured["rmse"] == pytest.approx(0.0, abs=1e-9) + assert measured["mean_observed"] == pytest.approx(measured["mean_simulated"]) + + @pytest.mark.unit + def test_the_efficiencies_are_the_ones_of_the_statistics(self, dataset): + # ``station_nse`` stays for the readers that only want the efficiency; + # it is taken from the statistics, so the two can never disagree. + evaluate(dataset.vector, dataset.context()) + + (record,) = dataset.records() + + assert record["station_nse"] == { + station: measured["nse"] for station, measured in record["station_metrics"].items() + } + + @pytest.mark.unit + def test_a_rejected_candidate_still_carries_the_moment_and_no_statistics(self, dataset): + from rubem.calibration.parameters import FREE_PARAMETERS + + vector = dataset.vector.copy() + vector[FREE_PARAMETERS.index("w_1")] = 0.6 + vector[FREE_PARAMETERS.index("w_2")] = 0.6 + + evaluate(vector, dataset.context()) + + (record,) = dataset.records() + assert record["error"] == "inadmissible" + assert record["station_metrics"] == {} + assert datetime.datetime.fromisoformat(record["started_at"]).tzinfo is not None + + @pytest.mark.unit + def test_a_station_outside_the_objective_is_measured_all_the_same(self, dataset): + # The objective averages station 1 alone; station 2 is reported anyway, + # which is the validation half of a calibration/validation split. + value = evaluate(dataset.vector, dataset.context(stations=("1",))) + + assert value == pytest.approx(0.0, abs=1e-9) + (record,) = dataset.records() + assert sorted(record["station_metrics"]) == ["1", "2"] + assert record["nse"] == pytest.approx(1.0, abs=1e-12) + + @pytest.mark.unit + def test_a_selection_that_names_no_shared_station_fails_the_evaluation(self, dataset): + from rubem.calibration.objective import INADMISSIBLE_OBJECTIVE + + value = evaluate(dataset.vector, dataset.context(stations=("nowhere",))) + + assert value == INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert "none of the stations" in record["error"] + + +class TestSimulateBest: + @pytest.mark.unit + def test_it_runs_the_candidate_and_keeps_the_table_it_wrote(self, dataset, tmp_path): + from rubem.calibration.parameters import vector_to_parameters + + output_dir = tmp_path / "best" + output_dir.mkdir() + + table = simulate_best( + dataset.document, + dataset.configuration.base_dir, + vector_to_parameters(dataset.vector), + str(output_dir), + "arn", + ) + + # The directory is left in place, which is the whole point: the series + # of the chosen candidate has to be readable after the run. + assert Path(table).is_file() + assert Path(table).parent == output_dir + assert sorted(entry.name for entry in output_dir.iterdir()) == [ + "metadata.json", + "tss_arn.csv", + ] + series = read_series(table) + assert sorted(series.stations) == ["1", "2"] + assert series.steps.tolist() == list(range(1, TIMESTEPS + 1)) + + @pytest.mark.unit + def test_it_reproduces_the_series_of_the_configuration_it_is_given(self, dataset, tmp_path): + from rubem.calibration.parameters import vector_to_parameters + + output_dir = tmp_path / "best" + output_dir.mkdir() + + table = simulate_best( + dataset.document, + dataset.configuration.base_dir, + vector_to_parameters(dataset.vector), + str(output_dir), + "arn", + ) + + simulated = read_series(table) + observed = read_series(dataset.observed_path) + for station, values in observed.stations.items(): + assert simulated.stations[station].tolist() == pytest.approx(values.tolist()) + + @pytest.mark.unit + def test_a_variable_the_run_does_not_write_raises(self, dataset, tmp_path): + from rubem.calibration.parameters import vector_to_parameters + + output_dir = tmp_path / "best" + output_dir.mkdir() + + # Unlike an evaluation, this one has nothing to rank a failure behind. + with pytest.raises(Exception, match="nope"): + simulate_best( + dataset.document, + dataset.configuration.base_dir, + vector_to_parameters(dataset.vector), + str(output_dir), + "nope", + ) diff --git a/tests/unit/test_cli_calibrate.py b/tests/unit/test_cli_calibrate.py new file mode 100644 index 00000000..984dcfd2 --- /dev/null +++ b/tests/unit/test_cli_calibrate.py @@ -0,0 +1,697 @@ +import json +from pathlib import Path + +import pytest + +from rubem.cli import main +from tests.helpers.synthetic import write_synthetic_dataset + + +@pytest.fixture(name="config_path") +def config_path_fixture(tmp_path): + path = tmp_path / "config.json" + path.write_text(json.dumps(write_synthetic_dataset(str(tmp_path))), encoding="utf8") + return path + + +@pytest.fixture(name="observed_path") +def observed_path_fixture(tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;1\n1;2.0\n2;3.0\n", encoding="utf8") + return path + + +class FakeCalibration: + """Stand in for ``rubem.calibration.runner.calibrate``. + + It records the arguments the command line built and answers with a result + whose files do not have to exist: the command only prints their paths. + """ + + def __init__(self, run_dir, best_nse=0.87): + from rubem.calibration.runner import CalibrationResult + + self.calls = [] + self.result = CalibrationResult( + best_parameters={ + "alpha": 4.5, + "beta": 0.5, + "w_1": 0.333, + "w_2": 0.333, + "w_3": 0.334, + "rcd": 5.0, + "f": 0.5, + "alpha_gw": 0.5, + "x": 0.5, + }, + best_nse=best_nse, + best_objective=169000.0, + evaluations=128, + generations=2, + success=True, + message="Optimization terminated successfully.", + run_dir=Path(run_dir), + evaluations_csv=Path(run_dir) / "evaluations.csv", + result_json=Path(run_dir) / "result.json", + calibrated_config=Path(run_dir) / "config-calibrated.json", + best_series=Path(run_dir) / "best_arn.csv", + stations_csv=Path(run_dir) / "stations.csv", + ) + + def __call__(self, config_path, observed_path, run_dir, settings): + self.calls.append((config_path, observed_path, run_dir, settings)) + return self.result + + +@pytest.fixture(name="fake_calibration") +def fake_calibration_fixture(tmp_path, monkeypatch): + fake = FakeCalibration(tmp_path / "calibration") + monkeypatch.setattr("rubem.calibration.runner.calibrate", fake) + return fake + + +class TestCliCalibrateArgumentsReachTheRunner: + @pytest.mark.unit + def test_the_paths_and_the_documented_defaults( + self, config_path, observed_path, tmp_path, fake_calibration, restore_logging + ): + run_dir = tmp_path / "calibration" + + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(run_dir), + ] + ) + + (given_config, given_observed, given_run_dir, settings) = fake_calibration.calls[0] + assert Path(given_config) == config_path + assert Path(given_observed) == observed_path + assert Path(given_run_dir) == run_dir + assert settings.variable == "arn" + assert settings.spinup_steps == 0 + assert settings.maxiter == 100 + assert settings.popsize == 15 + assert settings.seed is None + assert settings.temp_dir is None + assert settings.init == "sobol" + assert settings.strategy == "best1exp" + assert settings.polish is False + assert settings.bounds is None + assert settings.fixed is None + assert settings.stations is None + assert settings.allow_blocking_problems is False + + @pytest.mark.unit + def test_every_option_reaches_the_settings( + self, config_path, observed_path, tmp_path, fake_calibration, restore_logging + ): + temp_dir = tmp_path / "scratch" + + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + "--variable", + "rnf", + "--spinup-steps", + "12", + "--maxiter", + "3", + "--popsize", + "4", + "--seed", + "7", + "--workers", + "2", + "--temp-dir", + str(temp_dir), + ] + ) + + settings = fake_calibration.calls[0][3] + assert settings.variable == "rnf" + assert settings.spinup_steps == 12 + assert settings.maxiter == 3 + assert settings.popsize == 4 + assert settings.seed == 7 + assert settings.workers == 2 + assert Path(settings.temp_dir) == temp_dir + + @pytest.mark.unit + def test_the_runner_keeps_its_own_default_worker_count( + self, config_path, observed_path, tmp_path, fake_calibration, monkeypatch, restore_logging + ): + """Without ``--workers`` the command must not override the runner's default. + + The settings class is replaced by one whose default differs from the + ``None`` the option carries, so forwarding the option unconditionally + would be visible here. + """ + from dataclasses import dataclass + + from rubem.calibration.runner import CalibrationSettings + + @dataclass(frozen=True) + class SettingsWithAWorkerDefault(CalibrationSettings): + workers: int | None = 7 + + monkeypatch.setattr( + "rubem.calibration.runner.CalibrationSettings", SettingsWithAWorkerDefault + ) + + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + assert fake_calibration.calls[0][3].workers == 7 + + +class TestCliCalibrateOutput: + @pytest.mark.unit + def test_the_documented_progress_and_summary_are_printed( + self, config_path, observed_path, tmp_path, fake_calibration, capsys, restore_logging + ): + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + output = capsys.readouterr().out + result = fake_calibration.result + for expected in ( + "Loading configuration and validating inputs...", + "Calibration started...", + "Calibration finished successfully!", + "Best NSE: 0.87", + "Best objective: 169000", + "Evaluations: 128", + "Generations: 2", + str(result.evaluations_csv), + str(result.result_json), + str(result.calibrated_config), + f"Observed summary: {result.run_dir / 'observed.csv'}", + f"Stations table: {result.stations_csv}", + f"Best candidate series: {result.best_series}", + ): + assert expected in output, f"missing line: {expected!r}\n{output}" + + @pytest.mark.unit + def test_a_run_without_a_usable_nse_says_so( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + fake = FakeCalibration(tmp_path / "calibration", best_nse=None) + monkeypatch.setattr("rubem.calibration.runner.calibrate", fake) + + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + assert "Best NSE: n/a" in capsys.readouterr().out + + +class TestCliCalibrateFailures: + @pytest.mark.unit + def test_a_calibration_error_exits_with_one_without_a_traceback( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + from rubem.calibration.runner import CalibrationError + + def fail(*args, **kwargs): + raise CalibrationError("every evaluation failed: no such variable") + + monkeypatch.setattr("rubem.calibration.runner.calibrate", fail) + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + captured = capsys.readouterr().err + assert error.value.code == 1 + assert "every evaluation failed: no such variable" in captured + assert "Traceback" not in captured + + @pytest.mark.unit + def test_an_unexpected_failure_logs_one_traceback( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + def fail(*args, **kwargs): + raise RuntimeError("the pool died") + + monkeypatch.setattr("rubem.calibration.runner.calibrate", fail) + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + captured = capsys.readouterr().err + assert error.value.code == 1 + assert captured.count("Traceback (most recent call last)") == 1 + assert "RUBEM unexpectedly quit." in captured + + @pytest.mark.unit + def test_a_missing_scipy_exits_with_the_installation_guidance( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + monkeypatch.setattr("rubem._deps.missing_calibration_deps", lambda: ["scipy"]) + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + message = str(error.value) + assert "scipy" in message + assert 'pip install "rubem[calibration]"' in message + + @pytest.mark.unit + def test_an_interrupted_calibration_exits_with_two( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + def interrupt(*args, **kwargs): + raise KeyboardInterrupt + + monkeypatch.setattr("rubem.calibration.runner.calibrate", interrupt) + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + captured = capsys.readouterr().err + assert error.value.code == 2 + assert "RUBEM was interrupted by the user." in captured + assert "Traceback" not in captured + + @pytest.mark.unit + @pytest.mark.parametrize( + ("option", "value"), + [ + ("--maxiter", "0"), + ("--maxiter", "-1"), + ("--popsize", "0"), + ("--popsize", "-1"), + ("--workers", "0"), + ("--workers", "-1"), + ("--spinup-steps", "-1"), + ], + ) + def test_a_search_size_outside_its_range_is_rejected( + self, config_path, observed_path, tmp_path, fake_calibration, option, value, restore_logging + ): + """A generation count, a population and a worker count all start at one.""" + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + option, + value, + ] + ) + + assert error.value.code == 2 + assert not fake_calibration.calls + + @pytest.mark.unit + def test_a_negative_spin_up_is_rejected( + self, config_path, observed_path, tmp_path, fake_calibration, capsys, restore_logging + ): + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + "--spinup-steps", + "-1", + ] + ) + + assert error.value.code == 2 + assert not fake_calibration.calls + + @pytest.mark.unit + def test_a_missing_observed_file_is_rejected_before_anything_runs( + self, config_path, tmp_path, fake_calibration, capsys, restore_logging + ): + missing = tmp_path / "missing.csv" + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(missing), + "-o", + str(tmp_path / "calibration"), + ] + ) + + assert error.value.code == 2 + assert not fake_calibration.calls + assert "Loading configuration" not in capsys.readouterr().out + + @pytest.mark.unit + def test_an_incomplete_command_is_not_mapped_onto_the_legacy_spelling( + self, config_path, capsys, restore_logging + ): + """``calibrate`` is a subcommand, not the deprecated bare ``-c`` form.""" + with pytest.raises(SystemExit) as error: + main(["calibrate", "-c", str(config_path)]) + + captured = capsys.readouterr() + assert error.value.code == 2 + assert "Missing option '--observed'" in captured.err + assert "is deprecated" not in captured.err + + +class TestCliCalibrateHelp: + @pytest.mark.unit + def test_the_command_is_listed_in_the_help(self, capsys, restore_logging): + with pytest.raises(SystemExit) as error: + main(["--help"]) + + assert error.value.code == 0 + assert "calibrate" in capsys.readouterr().out + + @pytest.mark.unit + def test_the_help_documents_every_option(self, capsys, restore_logging): + with pytest.raises(SystemExit) as error: + main(["calibrate", "--help"]) + + assert error.value.code == 0 + output = capsys.readouterr().out + assert "Calibrate the model parameters against an observed series." in output + for option in ( + "--configfile", + "--observed", + "--run-dir", + "--variable", + "--spinup-steps", + "--maxiter", + "--popsize", + "--seed", + "--workers", + "--temp-dir", + "--bound", + "--fix", + "--stations", + "--init", + "--strategy", + "--polish", + "--no-polish", + ): + assert option in output, f"missing option: {option}\n{output}" + + @pytest.mark.unit + def test_importing_the_command_line_does_not_need_scipy(self): + """The CLI module stays importable without the optional extra.""" + import subprocess + import sys + + completed = subprocess.run( + [sys.executable, "-c", "import sys, rubem.cli; print('scipy' in sys.modules)"], + cwd=Path(__file__).resolve().parents[2], + capture_output=True, + text=True, + check=True, + ) + assert completed.stdout.strip() == "False" + + +class TestCliCalibrateDecisionSpaceOptions: + @pytest.mark.unit + def test_the_bounds_the_fixed_parameters_and_the_stations_reach_the_settings( + self, config_path, observed_path, tmp_path, fake_calibration, restore_logging + ): + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + "--bound", + "rcd=2:5", + "--bound", + "alpha=1:10", + "--fix", + "x=0", + "--stations", + "1, 2 ,3", + "--init", + "latinhypercube", + "--strategy", + "rand1bin", + "--polish", + ] + ) + + settings = fake_calibration.calls[0][3] + assert settings.bounds == {"rcd": (2.0, 5.0), "alpha": (1.0, 10.0)} + assert settings.fixed == {"x": 0.0} + # The ids are read as strings and the spacing around the commas is not + # part of them. + assert settings.stations == ("1", "2", "3") + assert settings.init == "latinhypercube" + assert settings.strategy == "rand1bin" + assert settings.polish is True + + @pytest.mark.unit + def test_no_polish_is_the_default_and_can_be_written_out( + self, config_path, observed_path, tmp_path, fake_calibration, restore_logging + ): + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + "--no-polish", + ] + ) + + assert fake_calibration.calls[0][3].polish is False + + @pytest.mark.unit + def test_allow_blocking_problems_reaches_the_settings( + self, config_path, observed_path, tmp_path, fake_calibration, restore_logging + ): + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + "--allow-blocking-problems", + ] + ) + + assert fake_calibration.calls[0][3].allow_blocking_problems is True + + @pytest.mark.unit + @pytest.mark.parametrize( + ("arguments", "message"), + [ + (["--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() diff --git a/tests/unit/test_deps.py b/tests/unit/test_deps.py index 410d5bca..41bcb5ae 100644 --- a/tests/unit/test_deps.py +++ b/tests/unit/test_deps.py @@ -18,3 +18,54 @@ def test_missing_dependencies_are_named_with_an_actionable_source(self, monkeypa assert name in message assert "conda-forge" in message assert _deps._ENVIRONMENT_YML_URL in message + + +class TestRequireCalibrationDeps: + @pytest.mark.unit + def test_present_dependencies_do_not_raise(self): + _deps.require_calibration_deps() + + @pytest.mark.unit + def test_nothing_is_missing_when_every_dependency_is_importable(self): + assert _deps.missing_calibration_deps() == [] + + @pytest.mark.unit + def test_an_unimportable_dependency_is_reported_as_missing(self, monkeypatch): + monkeypatch.setattr(_deps.importlib.util, "find_spec", lambda name: None) + + assert _deps.missing_calibration_deps() == list(_deps._CALIBRATION_DEPENDENCIES) + + @pytest.mark.unit + def test_a_broken_dependency_counts_as_missing(self, monkeypatch): + def raise_import_error(name): + raise ImportError(name) + + monkeypatch.setattr(_deps.importlib.util, "find_spec", raise_import_error) + + assert _deps.missing_calibration_deps() == list(_deps._CALIBRATION_DEPENDENCIES) + + @pytest.mark.unit + def test_the_message_names_the_dependency_and_the_extra(self): + message = _deps.calibration_deps_message(["scipy"]) + + assert "scipy" in message + assert 'pip install "rubem[calibration]"' in message + + @pytest.mark.unit + def test_the_message_defaults_to_every_calibration_dependency(self): + """The runner asks for the guidance without having computed the list.""" + assert _deps.calibration_deps_message() == _deps.calibration_deps_message( + list(_deps._CALIBRATION_DEPENDENCIES) + ) + + @pytest.mark.unit + def test_missing_dependencies_are_named_with_an_actionable_extra(self, monkeypatch): + monkeypatch.setattr(_deps.importlib.util, "find_spec", lambda name: None) + + with pytest.raises(SystemExit) as error: + _deps.require_calibration_deps() + + message = str(error.value) + for name in _deps._CALIBRATION_DEPENDENCIES: + assert name in message + assert 'pip install "rubem[calibration]"' in message