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 96344b8a..d432552f 100644 --- a/doc/source/api.rst +++ b/doc/source/api.rst @@ -122,6 +122,37 @@ Importing ``rubem.api`` does not require PCRaster or GDAL; running the model does. Without them, loading a configuration or running raises ``ImportError`` with the same installation guidance the command line prints. +Calibration +----------- + +``rubem calibrate`` is built on this API: each of its workers calls +``Model.run()`` in its own process, on a configuration derived from the one +being calibrated. The same search is available from Python as +``rubem.calibration.runner.calibrate``, which takes the configuration file, the +observed series, the run directory and a ``CalibrationSettings``: + +.. code-block:: python + + from rubem.calibration.runner import CalibrationSettings, calibrate + + if __name__ == "__main__": + result = calibrate( + "config.json", + "observed.csv", + "calibration", + CalibrationSettings(seed=42, maxiter=20, popsize=5), + ) + print(result.best_parameters, result.best_nse) + +The entry point guard is required: the calibration starts its workers with the +``spawn`` method, which imports the main module of the caller in every one of +them. SciPy is needed as well, from the ``rubem[calibration]`` extra. + +``rubem.calibration`` is **not** part of the stable surface yet: only +``rubem.api`` is, and the calibration package may change without a minor +version bump until it is promoted. The method behind the command is documented +in :doc:`Calibration `. + Logging ------- diff --git a/doc/source/calibration.rst b/doc/source/calibration.rst new file mode 100644 index 00000000..a8bb24cd --- /dev/null +++ b/doc/source/calibration.rst @@ -0,0 +1,555 @@ +Calibration +=========== + +.. role:: raw-html(raw) + :format: html + +The ``rubem calibrate`` command fits the calibration parameters of the model to +an observed streamflow series with a differential evolution search +[STORN1997]_. Every candidate the search proposes is a complete simulation: +the command builds the configuration of the candidate, runs it, samples the +resulting time series at the stations and scores it against the observations. + +The options of the command are listed in the user guide, under +:ref:`userguide:Running RUBEM`. This page documents what happens between them: +the objective, the parameters and their bounds, the settings of the search and +how to budget a run, the inputs it needs, what it writes, and how it uses the +machine it runs on. + +What the calibration fits +------------------------- + +The model has nine calibration parameters, the ``CALIBRATION`` section of the +configuration file: :math:`\alpha`, :math:`b`, :math:`w_1`, :math:`w_2`, +:math:`w_3`, :math:`RCD`, :math:`f`, :math:`\alpha_{GW}` and :math:`x`. They +are described in :ref:`overview:Calibration and Validation`. The search covers +eight of them; the slope factor weight :math:`w_3` is not searched but derived +from the other two weights, because the three weights of the potential runoff +coefficient must add up to 1. + +The observations are a streamflow series at the sample stations of the +configuration, one column per station. By default the simulated counterpart is +``arn``, the accumulated total runoff, in +:raw-html:`m3s-1`, which is the variable an observed +streamflow record is directly comparable with. ``--variable`` selects another +output variable, one of ``itp``, ``bfw``, ``srn``, ``eta``, ``lfw``, ``rec``, +``smc``, ``rnf`` and ``arn``; note that ``rnf`` is the total runoff of a single +cell, in millimeters, not a discharge, so an observed series compared with it +must be in the same quantity. + +Each evaluation runs the whole simulation period of the configuration. The +calibration never changes the period, the inputs or the aggregation: of what +the results depend on, it changes the nine parameters and nothing else. + +The objective function +---------------------- + +The efficiency +`````````````` + +The agreement between a simulated and an observed series at one station is the +Nash-Sutcliffe efficiency [MCCUEN2006]_: + +.. math:: + :label: calibration-nse + :nowrap: + + \[NSE = 1 - \frac{\sum_{t}{\left(Q_{sim,t} - Q_{obs,t}\right)^2}}{\sum_{t}{\left(Q_{obs,t} - \overline{Q_{obs}}\right)^2}}\] + +where: + +- :math:`Q_{sim,t}` – simulated value of the station at time step :math:`t`; +- :math:`Q_{obs,t}` – observed value of the station at time step :math:`t`; +- :math:`\overline{Q_{obs}}` – mean of the observed values used in the sum. + +Both sums run over the same pairs, and the mean is the mean of the observed +values of those pairs: a step dropped because the simulation has no value for +it is dropped from the mean as well, so numerator and denominator always +describe the same sample. + +Which values count +`````````````````` + +A pair counts only when both of its members count. A value does not count when +it is not finite (``NaN`` and the infinities), when it is ``-9999``, the gap +marker of the model's own tables and of the usual observation tables, or when +it is at or above ``1e30``. + +That last threshold is how the PCRaster missing value ``1e31`` is masked. The +value does not round-trip exactly through ``Float32`` — written and read back +it becomes ``9.999999...e30`` — so the mask rejects everything at or above +``1e30`` instead of comparing with the constant. No streamflow or runoff of a +real basin comes anywhere near that magnitude. An empty cell, or a cell that is +not a number at all, is read as ``NaN`` and masked the same way. + +Stations and the average +```````````````````````` + +The two series are aligned on the station ids they share and on the time steps +they share; a station or a step present on one side only is ignored. + +A station has no efficiency when fewer than two of its pairs count, or when all +of its valid observations are equal: a constant record has no variability +for the efficiency to explain and would divide by zero. Such a station is +reported as having no efficiency and is left out of the average, instead of +dragging it. The efficiency of a candidate is the mean over the stations that +do have one. + +A mismatch between the observations and the configuration is not a poor +candidate, and the two kinds of it that can be told before any model run are +refused before the search starts: an observed series with no time step in +common with the simulated steps that survive the spin-up, and one whose station +ids match none of the ids of the sample locations raster (an observed station +the configuration does not sample is reported once and ignored). With the +``zones`` aggregation the station ids are only known once a run has written +``zones_mapping.csv``, so that check is skipped. + +What cannot be told in advance surfaces inside the evaluations: a series whose +shared stations all lack an efficiency, because every observation is missing or +constant, fails each evaluation after its simulation has run. The evaluation is +then recorded with the message of the error in its ``error`` column and with +the objective ``1e30``, like any other failure, and the search carries on; the +calibration ends with an error only once no evaluation of the whole run has +succeeded. + +.. warning:: + + 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. + +The ninth parameter is derived rather than searched: + +.. math:: + :label: calibration-w3 + :nowrap: + + \[w_3 = 1 - \left(w_1 + w_2\right)\] + +so the search is eight-dimensional under the linear constraint +:math:`w_1 + w_2 \leq 1`. The constraint is handed to the optimizer, which +checks it before the objective and never spends a model run on a candidate +whose slope factor weight would come out negative. The derived weight must also +lie inside its own range, :math:`0 \leq w_3 \leq 1`; a candidate for which it +does not is rejected without a run and recorded with the error +``inadmissible``. Both checks are guards rather than a part of the search: the +optimizer keeps its trial vectors inside the bounds and drops the ones that +violate the constraint, so it does not propose an inadmissible candidate of its +own. + +The starting point of the search, ``x0``, is the parameter set of the +configuration being calibrated. The optimizer places it in the initial +population, so the configuration's own parameters are always among the +candidates of the first generation and a calibration can never return something +worse than what it started from without the table showing why. + +The search +---------- + +Settings +```````` + +The search is SciPy's ``differential_evolution`` with the following settings, +which the command does not expose because they describe the method rather than +the run: + +``strategy="best1exp"`` + The trial vector is built from the best member of the population plus a + scaled difference of two random members, with exponential crossover. A + best-based strategy exploits the best candidate found so far and converges + in fewer generations than the random-base ones [STORN1997]_, at some cost in + exploration, which the dithered mutation below offsets. Fewer generations is + what matters when one evaluation is a full simulation. + +``init="sobol"`` + The initial population is a Sobol' sequence over the bounds instead of a + random or Latin hypercube sample: it covers the eight-dimensional box far + more evenly, so the initial sample is not wasted on clusters and gaps. + +``mutation=(0.5, 1.0)`` + A pair, not a number: the mutation constant is drawn anew for every + generation between the two values (*dithering*). It widens the range of step + sizes the search takes without any tuning and reduces the chance of + stagnating in a local optimum. + +``recombination=0.7`` + The probability with which a trial vector keeps a component of the mutant + rather than of its parent. A high value suits parameters that interact, as + the three weights and the runoff coefficients do. + +``polish=False`` + No local refinement of the best candidate at the end. The local search would + spend model runs estimating gradients of an objective that is read back from + a simulated table, which does not provide them reliably. + +``updating="deferred"`` + The whole generation is evaluated before the population is updated. This is + what makes the generation a batch that can be spread over the worker + processes; the immediate alternative is inherently serial. + +``--seed`` seeds the random generator of the search: the same seed gives the +same initial population and the same random draws of the mutation and of the +crossover. That alone does not make two calibrations comparable candidate by +candidate, because from the first generation on the selection reads the +objectives, so the candidates a run proposes depend on the values the +simulations returned. Two runs follow the same path only when the simulation +itself is reproducible — see `A fixed drainage network`_. + +The evaluation budget +````````````````````` + +The number of population members is not ``--popsize``: it is + +.. math:: + :label: calibration-members + :nowrap: + + \[members = \max\left(5,\ popsize \times 8\right)\] + +with 8 the number of free parameters, rounded up to the next power of two, +because a Sobol' sequence is balanced only over a power-of-two sample. Every +generation evaluates one candidate per member, and the initial population is +one more round of evaluations on top of ``--maxiter`` generations, so + +.. math:: + :label: calibration-runs + :nowrap: + + \[runs \leq members \times \left(maxiter + 1\right)\] + +With the defaults, ``--popsize 15`` gives :math:`15 \times 8 = 120`, rounded up +to **128 members**, and ``--maxiter 100`` gives at most +:math:`128 \times 101 = 12\,928` model runs. A smaller search, ``--popsize 5``, +gives :math:`5 \times 8 = 40`, rounded up to **64 members**, so +``--popsize 5 --maxiter 20`` costs at most :math:`64 \times 21 = 1\,344` runs. + +Two things make the real count lower. Members whose weights violate +:math:`w_1 + w_2 \leq 1` are skipped by the optimizer without a model run and +leave no record. And the search stops as soon as the objective values of the +population agree: SciPy compares their standard deviation with its default +tolerance of 1% of their mean, and skips the comparison while any member of +the population is still infeasible. A search that converges this way stops +before ``--maxiter`` is reached. + +The number of members and the resulting upper bound are logged when the +calibration starts, so the budget can be checked before the machine is +committed to it: + +.. code-block:: console + + Calibrating 8 free parameter(s) with 128 population member(s) per generation and + at most 100 generation(s): up to 12928 model run(s), on 15 worker process(es). + +The size of the population the search actually built is recorded as +``population_size`` in :file:`result.json`, together with the number of +evaluations and of generations it really spent. + +Inputs +------ + +The observed series +``````````````````` + +``--observed`` takes one file, in either of two layouts, recognized from its +content. + +The **CSV layout** is the one the model writes for its own time series tables: +``;`` as the separator, a header whose first cell is ``0`` and whose remaining +cells are the station ids, and one row per time step carrying the step number +in the first column and one value per station: + +.. code-block:: console + + 0;1;2;3 + 1;12.4;3.1;0.8 + 2;15.9;4.0;1.2 + 3;-9999;4.4;1.0 + +The **PCRaster TSS layout** is accepted as well, both in its documented form (a +title line, the number of columns, one line per column name, then the rows) and +in the header-less form the model itself writes, whose stations are then +numbered ``1``, ``2``, ... in column order. + +.. note:: + + The station ids must be the ids of the sample locations of the + configuration. The easiest way to get the layout and the ids right is to run + the configuration once and use the table the run wrote as the template of + the observed file: :file:`tss_arn.csv` for the default ``point`` + aggregation, :file:`tss_arn_subcatchment.csv` or :file:`tss_arn_zones.csv` + for the other two. With ``zones`` the ids are the renumbered columns + ``1..N``, and the correspondence with the original zone ids is written to + :file:`zones_mapping.csv` in the output directory. + +The configuration must name the raster the aggregation samples from — the +sample locations raster for ``point`` and ``subcatchment``, the zones raster +for ``zones``. A configuration that does not name the raster its aggregation +needs is refused before any run: there would be no station to compare. + +A fixed drainage network +```````````````````````` + +When the configuration does not provide a Local Drain Direction raster, the +model derives one from the digital elevation model with PCRaster's +``lddcreate``. On the real basin DEMs that operation is not deterministic: two +runs of the same configuration produce drainage networks that differ in a few +cells, and the accumulated runoff ``arn``, which is routed over the network, +therefore differs between them. The other eight output variables are computed +cell by cell and are unaffected. + +A calibration on ``arn`` without a fixed network would compare series routed +over different networks from one evaluation to the next, and the differences +between candidates would be partly noise. **Calibrating on** ``arn`` +**requires** ``RASTERS.ldd`` **to name a fixed LDD raster.** Generate it once, +with the DEM of the basin, and keep it with the other inputs: + +.. code-block:: python + + import pcraster as pcr + + pcr.setclone("/basins/ipojuca/input/maps/clone.map") + dem = pcr.readmap("/basins/ipojuca/input/maps/dem.map") + ldd = pcr.lddcreate(dem, 1e31, 1e31, 1e31, 1e31) + pcr.report(ldd, "/basins/ipojuca/input/maps/ldd.map") + +.. code-block:: json + + { + "RASTERS": { + "ldd": "/basins/ipojuca/input/maps/ldd.map", + }, + } + +The same applies beyond ``arn`` when the time series are aggregated over +subcatchments: the subcatchment of each sample location is delineated over the +drainage network, so the areas the series are averaged over move with it, for +every variable. + +Artifacts of a run +------------------ + +The three artifacts of the run directory, ``-o``/``--run-dir``: + +:file:`evaluations.csv` + One row per evaluation, with the columns ``id``, ``pid``, ``alpha``, + ``beta``, ``w_1``, ``w_2``, ``w_3``, ``rcd``, ``f``, ``alpha_gw``, ``x``, + ``nse``, ``objective``, ``elapsed_seconds`` and ``error``. ``id`` identifies + the evaluation and ``pid`` the process that made it; ``nse`` is empty for a + candidate that produced no efficiency; ``error`` is empty for a successful + evaluation, ``inadmissible`` for a candidate rejected without a run, and the + type and message of the exception for a run that failed. The rows are + ordered by the candidate they evaluated, not by the moment they were + written, so that the table does not depend on how the parallel workers + happened to finish. + +:file:`result.json` + The summary: ``best_parameters`` (the nine parameters, the derived ``w_3`` + included), ``best_nse``, ``best_objective``, ``nfev`` and ``nit`` (the + evaluations and generations the search spent), ``success`` and ``message`` + from the optimizer, ``population_size``, and a ``settings`` object with + ``variable``, ``spinup_steps``, ``maxiter``, ``popsize``, ``seed``, + ``workers``, ``temp_dir``, ``init``, ``strategy``, ``mutation``, + ``recombination`` and ``polish``. The number of workers and the temporary + directory are the values the run resolved and used, not the defaults left + unset on the command line, so the summary describes a calibration that can + be repeated. + +:file:`-calibrated.json` + The calibrated configuration: the input document with the best parameters + in place, written in the format the input was written in, legacy or + format 1.0. Its paths are the ones the loader resolved, which are absolute + even when the input file wrote them relative to its own directory; moving + the file to another machine therefore means fixing the paths. + +Next to them, the :file:`evaluations/` directory holds one JSON record per +evaluation, written by the worker that made it. A record carries what the CSV +row carries and, in addition, ``station_nse``, the efficiency of each +individual station — which is where to look when the mean hides one station +behaving differently from the others. + +.. warning:: + + A run directory that already holds evaluation records of an earlier + calibration is refused before the search starts: consolidating two runs into + one table would mix them. Use an empty run directory for every calibration, + and keep the finished ones. + +Process model of the calibration +-------------------------------- + +The parent process loads the configuration once and validates its input files +once. The workers do not validate them again: they do not change during the +calibration, and re-reading every raster of the basin per evaluation would +dominate the cost. + +Every evaluation then runs in a **fresh interpreter**. PCRaster keeps its state +process-wide — ``setclone`` is global to the process and the raster memory of a +run is not reclaimed — so a worker reused for a second evaluation would carry +the state of the first. The pool is therefore started with the ``spawn`` method +and configured for one task per worker: a worker starts, evaluates one +candidate and exits. The generation is spread over as many such workers as +``--workers`` allows, which is also what makes the calibration parallel at all. + +Each worker builds the configuration of its candidate from the document of the +calibrated configuration: the nine parameters become the candidate's, the +output directory becomes a temporary directory of its own under ``--temp-dir``, +every raster series is disabled (an evaluation reads no raster of a previous +run, and writing them would dominate the cost of the run) and the time series +are reduced to the calibrated variable, as CSV. The aggregation the user +configured is kept, since it decides which areas the observed stations +correspond to. The temporary directory is removed when the evaluation ends, +whether it succeeded or not. + +The workers are silenced: a null handler on the root logger keeps the warnings +of a configuration that is not validated again out of the terminal, and +standard output is pointed at the null device so that the newline the PCRaster +framework writes when an interpreter exits is not printed once per evaluation. +Nothing is lost: a failure of an evaluation is recorded in that evaluation's +JSON record and in its row of the table, and the parent logs one summary at the +end saying how many evaluations failed and what the first failure was. One +failing candidate never ends a calibration; it is ranked behind every candidate +that was evaluated. + +Running on a cluster or a large machine +--------------------------------------- + +The cost of a calibration is the number of model runs times the cost of one +run, divided by the number of workers. Measure the second factor before +committing a machine: one ``rubem run`` of the configuration is an upper bound +on one evaluation, which writes no raster series and does not validate the +inputs again, and so costs a little less. On the Ipojuca basin, 228 monthly +time steps, an evaluation is of the order of half a minute. At 25 seconds +each, a default search on that basin, up to 12,928 runs on 15 workers, is on +the order of six hours; ``--popsize 5 --maxiter 20``, up +to 1,344 runs, is well under an hour. + +``--workers`` + Defaults to one less than the number of CPUs, which leaves a core for the + parent. On a shared machine, or under a batch scheduler, set it to the + number of cores the job was actually given rather than letting it read the + whole machine. Each worker holds one model run in memory, plus the parent, + so the memory of the job is roughly ``--workers`` times the memory of a + single ``rubem run`` of the same configuration; that, and not the core + count, is usually what limits the number of workers on a large basin. + +``--temp-dir`` + The per-evaluation output directories are created here. Point it at fast + local storage of the node — the scratch directory of the job, not a shared + network filesystem: every evaluation creates a directory, writes a small + table into it and removes it again, and a network filesystem turns that into + the slowest part of the evaluation. The run directory itself may stay on + shared storage; it receives one small JSON record per evaluation. + +Shorten the simulation period first + An exploratory calibration over a few months of ``SIM_TIME``, with a small + ``--popsize`` and ``--maxiter``, costs a fraction of the full search and + still shows whether the observed series, the station ids and the spin-up are + set up as intended. Move to the full period only once a short run finishes + and produces a sensible efficiency. + +Continue from a calibrated configuration + :file:`-calibrated.json` is an ordinary configuration file. Passing + it back to ``rubem calibrate`` as ``-c`` starts a second search from the + parameters the first one found, since the starting point of a search is the + parameter set of the configuration it is given — useful to refine a result + with a different seed, a longer period or a larger population. + +Keep the run directory + :file:`evaluations.csv` and :file:`result.json` are the record of how a + parameter set was obtained: the settings, the seed, the number of + evaluations and the whole response surface the search sampled. They are + small, and they are what makes a published calibration reproducible. diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index e39985a1..da34c43b 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -13,6 +13,23 @@ The format follows `Keep a Changelog `__. Added ````` +- Added the calibration page of the documentation (the objective, the free + parameters and their bounds, the derivation of ``w3``, the differential + evolution settings and the evaluation budget, the observed series formats, + the fixed drainage network a calibration on ``arn`` needs, the artifacts, the + process model and hints for large machines), the ``dataset`` pytest marker + with the ``RUBEM_DATASET_DIR`` convention for tests that need the published + datasets, and a calibration smoke test on the Ipojuca basin + (`#345 `__). +- Added ``rubem calibrate``: a differential evolution (SciPy) over the eight free + calibration parameters, ``w3`` derived from the other two weights, minimizing + ``1000 (100 (1 - NSE))^2`` on the station series of one output variable + (``arn`` by default) against an observed CSV or PCRaster TSS file, with + masked values, a zero-variance guard and a spin-up window; every candidate + runs in a fresh spawned worker through the Python API, records one JSON + file, and the run directory receives ``evaluations.csv``, ``result.json`` + and ``-calibrated.json``. SciPy comes with the optional extra + ``rubem[calibration]`` (`#343 `__). - Added ``rubem.api``, the public Python surface: ``Model.from_file`` and ``Model.from_config`` load a configuration file or document, ``Model.run()`` runs the simulation in the current process and ``Model.run_isolated()`` in a 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. + -h, --help Show this message and exit. + +``calibrate`` needs SciPy, which comes with ``pip install +"rubem[calibration]"``. It writes ``evaluations.csv``, ``result.json`` +and ``-calibrated.json`` into the ``--run-dir``. The method itself (the +objective, the parameter bounds, the derivation of ``w3``, how to budget the +number of evaluations and what a run needs from the configuration) is +documented in :doc:`Calibration `. + +.. code-block:: console + + $ rubem calibrate -c project-config.json --observed observed.csv -o calibration + Loading configuration and validating inputs... + Calibration started... + Calibration finished successfully! + Best NSE: 0.873210 + Best objective: 160757 + Evaluations: 2048 + Generations: 16 + Evaluations table: /basins/ipojuca/calibration/evaluations.csv + Result: /basins/ipojuca/calibration/result.json + Calibrated configuration: /basins/ipojuca/calibration/project-config-calibrated.json + Use ``-V`` or ``--version`` to get the version of the RUBEM. .. code-block:: console diff --git a/pyproject.toml b/pyproject.toml index 14df2be4..29b24f6a 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", @@ -90,6 +95,7 @@ markers = [ "unit: mark a test as a unit test.", "integration: mark a test as a integration test.", "exact: mark a test that reproduces the goldens byte-for-byte on the golden environment.", + "dataset: mark a test that needs the published datasets (RUBEM_DATASET_DIR).", ] [tool.coverage.run] diff --git a/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..edf32e89 --- /dev/null +++ b/rubem/calibration/_worker.py @@ -0,0 +1,281 @@ +"""Evaluation of one candidate, the entry point the calibration workers run. + +The module lives at the top level of the package and keeps its entry point, +:func:`evaluate`, at module level, because the process pool of the calibration +starts its workers with the ``spawn`` method: the child imports this module by +name to find the function it has to call, so the function must be picklable by +reference and everything it needs must travel as plain data in an +:class:`EvaluationContext`. + +One evaluation builds the configuration of the candidate from the document of +the calibrated configuration, runs the model in the current process, reads the +time series it wrote and compares it with the observed series. It never calls +:meth:`rubem.api.Model.run_isolated`: the worker is already the isolation layer, +and an isolated run inside it would nest one process pool in another. + +Nothing here raises: a candidate outside the admissible region is rejected +without a run, and any failure of the run or of the evaluation is recorded and +turned into :data:`rubem.calibration.objective.INADMISSIBLE_OBJECTIVE`, so that +one failing candidate does not end the calibration. +""" + +import copy +import json +import logging +import os +import shutil +import tempfile +import time +import uuid +from collections.abc import Sequence +from dataclasses import dataclass +from pathlib import Path + +import numpy as np + +from ..configuration.model_configuration_file_v1 import VARIABLE_IDS +from .objective import ( + INADMISSIBLE_OBJECTIVE, + Series, + evaluate_series, + objective, + read_series, +) +from .parameters import is_admissible, vector_to_parameters + +logger = logging.getLogger(__name__) + +_INADMISSIBLE_ERROR = "inadmissible" +"""The error recorded for a candidate that was rejected without being run.""" + + +@dataclass(frozen=True) +class EvaluationContext: + """Everything an evaluation needs, as data that crosses to a worker process. + + :param document: The calibrated configuration as a format 1.0 document, + with its paths already anchored. + :type document: dict + + :param base_dir: Directory the relative paths of the document are anchored + on, ``None`` when they are absolute already. + :type base_dir: str | None + + :param variable: Id of the output variable the efficiency is computed on + (``arn``, the accumulated total runoff, by default). + :type variable: str + + :param observed: The observed series, read once by the parent. + :type observed: rubem.calibration.objective.Series + + :param spinup_steps: Number of initial time steps excluded from the + efficiency. + :type spinup_steps: int + + :param temp_dir: Directory the per-evaluation output directories are + created in and removed from. + :type temp_dir: str + + :param evaluations_dir: Directory the JSON record of every evaluation is + written to. + :type evaluations_dir: str + + :param validate_input: Whether every evaluation revalidates the input files + and their content. ``False`` by default: the parent validated them once + before the calibration started, and they do not change during it. + :type validate_input: bool + """ + + document: dict + base_dir: str | None + variable: str + observed: Series + spinup_steps: int + temp_dir: str + evaluations_dir: str + validate_input: bool = False + + +def evaluate(vector: Sequence[float] | np.ndarray, context: EvaluationContext) -> float: + """Run one candidate and return the objective value of its simulation. + + A candidate outside the admissible region is rejected before any model run + and recorded with the error ``inadmissible``. Otherwise the model is run in + the current process, into a temporary output directory that is removed + again, and the efficiency of the resulting time series is turned into the + objective. Every evaluation, successful or not, leaves one JSON record in + the evaluations directory of the context. + + :param vector: The candidate, the eight free values in the order of + :data:`rubem.calibration.parameters.FREE_PARAMETERS`. + :type vector: collections.abc.Sequence[float] | numpy.ndarray + + :param context: What the evaluation runs against. + :type context: EvaluationContext + + :return: The objective value, or + :data:`rubem.calibration.objective.INADMISSIBLE_OBJECTIVE` when the + candidate was rejected or its evaluation failed. + :rtype: float + """ + started = time.perf_counter() + try: + parameters = vector_to_parameters(vector) + except ValueError as error: + _write_record(context, {}, None, {}, INADMISSIBLE_OBJECTIVE, 0.0, _describe(error)) + return INADMISSIBLE_OBJECTIVE + + if not is_admissible(vector): + logger.debug("Rejecting an inadmissible candidate without running the model.") + _write_record( + context, + parameters, + None, + {}, + INADMISSIBLE_OBJECTIVE, + time.perf_counter() - started, + _INADMISSIBLE_ERROR, + ) + return INADMISSIBLE_OBJECTIVE + + output_dir = tempfile.mkdtemp(prefix="evaluation-", dir=context.temp_dir) + try: + nse, station_nse = _run_and_evaluate(context, parameters, output_dir) + except Exception as error: + # One failing candidate must not end the calibration: it is recorded and + # ranked behind every candidate that was actually evaluated. + logger.warning("The evaluation of a candidate failed: %s", error) + _write_record( + context, + parameters, + None, + {}, + INADMISSIBLE_OBJECTIVE, + time.perf_counter() - started, + _describe(error), + ) + return INADMISSIBLE_OBJECTIVE + finally: + shutil.rmtree(output_dir, ignore_errors=True) + + value = objective(nse) + _write_record(context, parameters, nse, station_nse, value, time.perf_counter() - started, None) + return value + + +def _run_and_evaluate( + context: EvaluationContext, + parameters: dict[str, float], + output_dir: str, +) -> tuple[float, dict[str, float | None]]: + """Run the candidate and return its mean efficiency and the efficiency per station.""" + from ..api import Model + + document = _derived_document(context, parameters, output_dir) + result = Model.from_config( + document, validate_input=context.validate_input, base_dir=context.base_dir + ).run() + simulated = read_series(_time_series_table(result, context.variable)) + return evaluate_series(simulated, context.observed, context.spinup_steps) + + +def _time_series_table(result, variable: str) -> Path: + """Return the CSV table the run wrote for ``variable``. + + :raises RuntimeError: If the run reports no table for the variable. + """ + tables = result.time_series.get(variable, ()) + for path in tables: + if path.suffix.lower() == ".csv": + return path + raise RuntimeError( + f"The run wrote no CSV time series for the variable '{variable}'; it reports " + f"{', '.join(str(path) for path in tables) or 'no table'}." + ) + + +def _derived_document( + context: EvaluationContext, + parameters: dict[str, float], + output_dir: str, +) -> dict: + """Return the configuration document of one candidate. + + The calibrated document is copied and four things are changed: the + calibration parameters become the candidate's, the outputs go to the + temporary directory of the evaluation, every raster series is disabled (an + evaluation reads no raster, and writing them would dominate its cost) and + the time series are reduced to the calibrated variable, as CSV. The + aggregation the user configured is kept, since it decides which areas the + observed stations correspond to. + """ + document = copy.deepcopy(context.document) + document["model_calibration_parameters"] = { + "alpha": parameters["alpha"], + "b": parameters["beta"], + "w_1": parameters["w_1"], + "w_2": parameters["w_2"], + "w_3": parameters["w_3"], + "rcd": parameters["rcd"], + "f": parameters["f"], + "alpha_gw": parameters["alpha_gw"], + "x": parameters["x"], + } + output = document.setdefault("model_simulation_output", {}) + previous_rasters = output.get("raster_series", {}) + previous_samples = output.get("time_series_samples", {}) + output["dir_path"] = output_dir + output["raster_series"] = { + "no_data_value": previous_rasters.get("no_data_value", -9999), + "formats": [], + } + output["time_series_samples"] = { + **dict.fromkeys(VARIABLE_IDS, False), + context.variable: True, + "formats": ["CSV"], + "aggregation": previous_samples.get("aggregation", "point"), + } + return document + + +def _describe(error: BaseException) -> str: + """Return the recorded description of an exception: its type and its message.""" + return f"{type(error).__name__}: {error}" + + +def _write_record( + context: EvaluationContext, + parameters: dict[str, float], + nse: float | None, + station_nse: dict[str, float | None], + value: float, + elapsed_seconds: float, + error: str | None, +) -> None: + """Write the JSON record of one evaluation, atomically. + + Every worker writes into the same directory while the parent may be reading + it, so a record is written to a temporary file next to its destination and + renamed onto it: a reader either does not see the record yet or sees it + whole. The name is a fresh UUID, so two workers never collide. + """ + record = { + "id": uuid.uuid4().hex, + "pid": os.getpid(), + "parameters": parameters, + "nse": nse, + "station_nse": station_nse, + "objective": value, + "elapsed_seconds": elapsed_seconds, + "error": error, + } + directory = Path(context.evaluations_dir) + target = directory / f"{record['id']}.json" + handle, temporary = tempfile.mkstemp(prefix=f".{target.name}.", suffix=".tmp", dir=directory) + try: + with os.fdopen(handle, "w", encoding="utf-8") as file: + json.dump(record, file) + Path(temporary).replace(target) + except BaseException: + Path(temporary).unlink(missing_ok=True) + raise diff --git a/rubem/calibration/objective.py b/rubem/calibration/objective.py new file mode 100644 index 00000000..7839851e --- /dev/null +++ b/rubem/calibration/objective.py @@ -0,0 +1,372 @@ +"""The calibration objective: Nash-Sutcliffe efficiency over the sample stations. + +The module reads the observed and the simulated time series, masks the values +that stand for a missing observation, computes the Nash-Sutcliffe efficiency +(NSE) of each station and turns the mean efficiency into the objective the +optimizer minimizes. + +Nothing here imports SciPy: the efficiency and the objective are evaluated with +NumPy alone. +""" + +import csv +import logging +from dataclasses import dataclass +from pathlib import Path + +import numpy as np + +from .._paths import PathInput, as_path + +logger = logging.getLogger(__name__) + +MISSING_VALUES = (-9999.0,) +"""Values written by the model and by the usual observation tables for a gap.""" + +PCRASTER_MISSING = 1e31 +"""The PCRaster missing value. + +Written as ``Float32`` it does not round-trip exactly (``float(np.float32(1e31))`` +is ``9.999999...e30``), so the mask rejects everything at or above ``1e30`` +instead of comparing with this constant: no streamflow or runoff of a real +basin comes anywhere near that magnitude. +""" + +_MISSING_THRESHOLD = 1e30 + +INADMISSIBLE_OBJECTIVE = 1.0e30 +"""The objective of a candidate that is not run. + +The worst objective a run can produce is bounded by the efficiencies it can +reach: an NSE of -100, already an absurd simulation, gives +``1000 * (100 * 101) ** 2``, about ``1.0e11``. ``1.0e30`` is far above any such +value, so an inadmissible or failed candidate is always ranked behind every +candidate that was actually evaluated. It is finite so that it sorts, is +written to JSON and lands in the table of evaluations like any other objective. +""" + + +def valid_mask(values: np.ndarray) -> np.ndarray: + """Return the mask of the values that count as an observation. + + A value counts when it is finite, is not one of :data:`MISSING_VALUES` and + is below ``1e30``, which covers the PCRaster missing value in any of its + round-trips. + + :param values: Values of one station, in step order. + :type values: numpy.ndarray + + :return: A boolean array of the shape of ``values``. + :rtype: numpy.ndarray + """ + array = np.asarray(values, dtype=np.float64) + mask = np.isfinite(array) & (array < _MISSING_THRESHOLD) + for missing in MISSING_VALUES: + mask &= array != missing + return mask + + +def nash_sutcliffe(simulated: np.ndarray, observed: np.ndarray) -> float | None: + """Return the Nash-Sutcliffe efficiency of one station, or ``None``. + + The efficiency is computed over the pairs whose two members are valid, from + its definition:: + + NSE = 1 - sum((simulated - observed) ** 2) / sum((observed - mean(observed)) ** 2) + + Both sums run over the same pairs, and the mean is the mean of the observed + values of those pairs: a step dropped because the simulation has no value + for it is dropped from the mean as well, so that the numerator and the + denominator always describe the same sample. + + :param simulated: Simulated values, in step order. + :type simulated: numpy.ndarray + + :param observed: Observed values, aligned with ``simulated``. + :type observed: numpy.ndarray + + :return: The efficiency, or ``None`` when fewer than two pairs are valid or + when the valid observed values are all equal, which leaves the + denominator at zero and the efficiency undefined. + :rtype: float | None + + :raises ValueError: If the two series do not have the same length. + """ + simulated = np.asarray(simulated, dtype=np.float64) + observed = np.asarray(observed, dtype=np.float64) + if simulated.shape != observed.shape: + raise ValueError( + "The simulated and the observed series must have the same shape, got " + f"{simulated.shape} and {observed.shape}." + ) + + mask = valid_mask(simulated) & valid_mask(observed) + if int(np.count_nonzero(mask)) < 2: + return None + + simulated_valid = simulated[mask] + observed_valid = observed[mask] + # The zero-variance guard: a station whose observations are all equal has no + # variability for the efficiency to explain. Constancy is tested on the + # values themselves, since the sum of the squared deviations of a constant + # series is only approximately zero in floating point. + if np.all(observed_valid == observed_valid[0]): + return None + + residuals = float(np.sum((simulated_valid - observed_valid) ** 2)) + variability = float(np.sum((observed_valid - np.mean(observed_valid)) ** 2)) + return 1.0 - residuals / variability + + +def objective(nse: float) -> float: + """Return the objective value of an efficiency, ``1000 (100 (1 - NSE))^2``. + + The objective is zero for a perfect simulation and grows quadratically as + the efficiency falls, which is what the differential evolution minimizes. + + :param nse: The Nash-Sutcliffe efficiency of the candidate. + :type nse: float + + :return: The objective value. + :rtype: float + """ + return 1000.0 * (100.0 * (1.0 - nse)) ** 2 + + +@dataclass(frozen=True) +class Series: + """A time series table: the steps and the values of each station. + + :param steps: The time steps, one entry per row, as integers. + :type steps: numpy.ndarray + + :param stations: Station id to the values of that station, as ``float64`` + arrays aligned with :attr:`steps`. + :type stations: dict[str, numpy.ndarray] + """ + + steps: np.ndarray + stations: dict[str, np.ndarray] + + +def _as_value(token: str) -> float: + """Return ``token`` as a float, or ``NaN`` when it is empty or not a number. + + An entry a table leaves empty or writes in a form that is not a number is a + gap; it becomes ``NaN`` and the mask drops it, as it drops the gaps written + as ``-9999`` or as the PCRaster missing value. + """ + try: + return float(token) + except ValueError: + return float("nan") + + +def _as_step(token: str, path: Path, number: int) -> int: + """Return the time step written in the first column of a row. + + :raises ValueError: If the token is not a number, naming the file and the + row, since a table whose first column is a date or a label is a common + mistake and the bare conversion error does not say where it came from. + """ + try: + return int(float(token)) + except ValueError: + raise ValueError( + f"Row {number} of the time series file {path} does not start with a time " + f"step: {token!r} is not a number." + ) from None + + +def _build(ids: list[str], steps: list[int], rows: list[list[str]], path: Path) -> Series: + """Assemble a :class:`Series` from the parsed header, steps and value rows.""" + duplicated = sorted({station for station in ids if ids.count(station) > 1}) + if duplicated: + raise ValueError( + f"The time series file {path} has repeated station ids: {', '.join(duplicated)}." + ) + if not ids: + raise ValueError(f"The time series file {path} has no station column.") + if not rows: + raise ValueError(f"The time series file {path} has no data rows.") + + values = np.array( + [[_as_value(token) for token in row] for row in rows], + dtype=np.float64, + ) + return Series( + steps=np.asarray(steps, dtype=np.int64), + stations={station: values[:, column] for column, station in enumerate(ids)}, + ) + + +def _read_csv(path: Path) -> Series: + """Read the CSV layout of the model's time series tables. + + The separator is ``;``, the header is ``0;;...`` and every row + carries the time step in the first column and one value per station. + """ + with path.open(encoding="utf8", newline="") as table: + records = [ + row for row in csv.reader(table, delimiter=";") if any(cell.strip() for cell in row) + ] + if not records: + raise ValueError(f"The time series file {path} is empty.") + + header = [cell.strip() for cell in records[0]] + ids = header[1:] + steps: list[int] = [] + rows: list[list[str]] = [] + for number, record in enumerate(records[1:], 2): + cells = [cell.strip() for cell in record] + if len(cells) != len(header): + raise ValueError( + f"Row {number} of the time series file {path} has {len(cells)} column(s), " + f"the header gives {len(header)}." + ) + steps.append(_as_step(cells[0], path, number)) + rows.append(cells[1:]) + return _build(ids, steps, rows, path) + + +def _read_tss(path: Path) -> Series: + """Read a PCRaster time series (``.tss``) file. + + The layout is a title line, a line with the number of columns, one line per + column name starting with the time step column, and then one whitespace + separated row per time step. The model writes its own ``.tss`` files + without that header (see ``rubem.file._file_conversions.tss2csv``, which + reads every line of them as data); such a file is read as well, and its + stations are then numbered ``1``, ``2``, ... in column order, which is the + numbering of the sample locations of a usual dataset. + """ + lines = [line for line in path.read_text(encoding="utf8").splitlines() if line.strip()] + if not lines: + raise ValueError(f"The time series file {path} is empty.") + + ids: list[str] | None = None + data_start = 0 + if len(lines) > 1 and lines[1].strip().isdigit(): + columns = int(lines[1].strip()) + names = [line.strip() for line in lines[2 : 2 + columns]] + if len(names) == columns and columns >= 2: + ids = names[1:] + data_start = 2 + columns + + rows: list[list[str]] = [] + steps: list[int] = [] + for number, line in enumerate(lines[data_start:], data_start + 1): + cells = line.split() + if ids is None: + ids = [str(column) for column in range(1, len(cells))] + if len(cells) != len(ids) + 1: + raise ValueError( + f"Row {number} of the time series file {path} has {len(cells)} column(s), " + f"the column names give {len(ids) + 1}." + ) + steps.append(_as_step(cells[0], path, number)) + rows.append(cells[1:]) + return _build(ids or [], steps, rows, path) + + +def read_series(path: PathInput) -> Series: + """Read a time series table, in the model's CSV layout or as a PCRaster TSS. + + The layout is recognised from the content: a first line carrying the ``;`` + separator is the CSV table the model writes, anything else is read as a + PCRaster time series. Missing entries are kept as they are written; the + mask of :func:`valid_mask` is what decides whether a value counts. + + :param path: The file to read. + :type path: str | os.PathLike[str] + + :return: The parsed series. + :rtype: Series + + :raises FileNotFoundError: If the file does not exist. + :raises ValueError: If the file is empty, has no station column, repeats a + station id or has a row that does not match its header. + """ + file_path = as_path(path) + text = file_path.read_text(encoding="utf8") + first_line = next((line for line in text.splitlines() if line.strip()), "") + if ";" in first_line: + return _read_csv(file_path) + return _read_tss(file_path) + + +def evaluate_series( + simulated: Series, + observed: Series, + spinup_steps: int = 0, +) -> tuple[float, dict[str, float | None]]: + """Return the mean efficiency of a run and the efficiency of each station. + + The two series are aligned on the steps they share, after the steps up to + ``spinup_steps`` have been dropped from both, and on the station ids they + share: a station or a step present on one side only is ignored. The mean is + taken over the stations that have an efficiency, so that a station of zero + variance, reported as ``None``, does not drag the mean. + + A step repeated in a table is aligned on its first row, since that is the + occurrence :func:`numpy.intersect1d` reports for a repeated value. + + :param simulated: The series the run wrote. + :type simulated: Series + + :param observed: The observed series. + :type observed: Series + + :param spinup_steps: Number of initial steps to exclude, defaults to ``0``. + :type spinup_steps: int, optional + + :return: The mean efficiency and the efficiency of every shared station. + :rtype: tuple[float, dict[str, float | None]] + + :raises ValueError: If the two series share no station id, share no time + step, or if no shared station has an efficiency. + """ + shared_ids = sorted(set(simulated.stations) & set(observed.stations)) + if not shared_ids: + raise ValueError( + "The simulated and the observed series have no station in common: the " + f"simulation has {', '.join(sorted(simulated.stations)) or 'none'} and the " + f"observations have {', '.join(sorted(observed.stations)) or 'none'}." + ) + + simulated_steps = np.asarray(simulated.steps, dtype=np.int64) + observed_steps = np.asarray(observed.steps, dtype=np.int64) + simulated_kept = simulated_steps > spinup_steps + observed_kept = observed_steps > spinup_steps + shared_steps, simulated_rows, observed_rows = np.intersect1d( + simulated_steps[simulated_kept], + observed_steps[observed_kept], + return_indices=True, + ) + if shared_steps.size == 0: + raise ValueError( + "The simulated and the observed series have no time step in common after " + f"the {spinup_steps} spin-up step(s): the simulation has " + f"{simulated_steps[simulated_kept].tolist()} and the observations have " + f"{observed_steps[observed_kept].tolist()}." + ) + + simulated_index = np.flatnonzero(simulated_kept)[simulated_rows] + observed_index = np.flatnonzero(observed_kept)[observed_rows] + + station_nse: dict[str, float | None] = { + station: nash_sutcliffe( + simulated.stations[station][simulated_index], + observed.stations[station][observed_index], + ) + for station in shared_ids + } + + values = [nse for nse in station_nse.values() if nse is not None] + if not values: + raise ValueError( + "No station has a Nash-Sutcliffe efficiency: every shared station " + f"({', '.join(shared_ids)}) has fewer than two valid pairs or observations " + "without variance." + ) + return float(np.mean(values)), station_nse diff --git a/rubem/calibration/parameters.py b/rubem/calibration/parameters.py new file mode 100644 index 00000000..de14d6ac --- /dev/null +++ b/rubem/calibration/parameters.py @@ -0,0 +1,184 @@ +"""The decision vector of the calibration and its translation to model parameters. + +The differential evolution searches eight free parameters, ``alpha``, ``beta``, +``w_1``, ``w_2``, ``rcd``, ``f``, ``alpha_gw`` and ``x``; the slope factor +weight ``w_3`` is not searched but derived from the other two weights, since the +three weights of the weighted runoff coefficient must add up to 1. The search is +therefore eight-dimensional with the linear constraint ``w_1 + w_2 <= 1``, which +keeps the derived weight inside its own range. + +The bounds are the ones the application settings declare (``appsettings.json``, +section ``value_ranges.variables``), the same ranges +:class:`rubem.configuration.calibration_parameters.CalibrationParameters` +validates against, so a candidate the optimizer proposes inside the bounds is a +configuration the model accepts. + +Nothing here imports SciPy at module level; :func:`weights_constraint`, the only +function that needs it, imports it when it is called. +""" + +import logging +from collections.abc import Mapping, Sequence +from typing import TYPE_CHECKING + +import numpy as np + +from ..configuration._ranges import variable_range + +if TYPE_CHECKING: + from scipy.optimize import LinearConstraint + +logger = logging.getLogger(__name__) + +FREE_PARAMETERS = ("alpha", "beta", "w_1", "w_2", "rcd", "f", "alpha_gw", "x") +"""The parameters the optimizer searches, in the order of the decision vector.""" + +DERIVED_PARAMETER = "w_3" +"""The slope factor weight, derived from ``w_1`` and ``w_2`` instead of searched.""" + +CALIBRATION_PARAMETERS = ( + "alpha", + "beta", + "w_1", + "w_2", + "w_3", + "rcd", + "f", + "alpha_gw", + "x", +) +"""The nine calibration parameters of the model, in the order of the configuration.""" + +_ALIASES = {"b": "beta", "w1": "w_1", "w2": "w_2", "w3": "w_3"} +"""Spellings the legacy configuration file accepts for a calibration parameter.""" + + +def bounds() -> list[tuple[float, float]]: + """Return the ``(minimum, maximum)`` bound of every free parameter. + + The bounds come from the application settings, in the order of + :data:`FREE_PARAMETERS`, which is the order of the decision vector. + + :return: One ``(minimum, maximum)`` pair per free parameter. + :rtype: list[tuple[float, float]] + """ + return [variable_range(name) for name in FREE_PARAMETERS] + + +def vector_to_parameters(vector: Sequence[float] | np.ndarray) -> dict[str, float]: + """Return the nine calibration parameters a decision vector stands for. + + The eight free values are read in the order of :data:`FREE_PARAMETERS` and + the slope factor weight is derived as ``w_3 = 1 - w_1 - w_2``, so that the + three weights always add up to 1. + + :param vector: The eight free values, in the order of :data:`FREE_PARAMETERS`. + :type vector: collections.abc.Sequence[float] | numpy.ndarray + + :return: The nine parameters, keyed by the names of + :class:`rubem.configuration.calibration_parameters.CalibrationParameters`. + :rtype: dict[str, float] + + :raises ValueError: If the vector does not have eight entries. + """ + values = np.asarray(vector, dtype=np.float64).ravel() + if values.size != len(FREE_PARAMETERS): + raise ValueError( + f"The decision vector must have {len(FREE_PARAMETERS)} entries " + f"({', '.join(FREE_PARAMETERS)}), got {values.size}." + ) + parameters = {name: float(value) for name, value in zip(FREE_PARAMETERS, values, strict=True)} + # The two weights are added first and the sum is subtracted once. Taking + # them away one at a time, ``1 - w_1 - w_2``, rounds twice: for 20 of the + # 101 pairs of two-decimal weights that add up to 1, among them + # ``w_1 = 0.33, w_2 = 0.67``, it lands one ulp below zero, and + # :class:`rubem.configuration.calibration_parameters.CalibrationParameters` + # then refuses the run. With a single subtraction the weight is exactly + # zero at the boundary and never negative while ``w_1 + w_2 <= 1``. + parameters[DERIVED_PARAMETER] = 1.0 - (parameters["w_1"] + parameters["w_2"]) + return {name: parameters[name] for name in CALIBRATION_PARAMETERS} + + +def parameters_to_vector(parameters: Mapping[str, float]) -> np.ndarray: + """Return the decision vector of a set of calibration parameters. + + This is how the configuration under calibration becomes the ``x0`` of the + differential evolution. Both spellings the configuration files use are + accepted, the canonical ``beta``, ``w_1``, ``w_2`` and the legacy ``b``, + ``w1``, ``w2``; the derived ``w_3`` is ignored, since it is not searched. + The straightforward way to obtain the mapping is to dump the loaded + :class:`rubem.configuration.calibration_parameters.CalibrationParameters` + of a configuration, which has already normalised the spellings. + + :param parameters: The calibration parameters, by name. + :type parameters: collections.abc.Mapping[str, float] + + :return: The eight free values, in the order of :data:`FREE_PARAMETERS`. + :rtype: numpy.ndarray + + :raises KeyError: If a free parameter is missing from the mapping. + """ + normalised = {_ALIASES.get(name, name): value for name, value in parameters.items()} + missing = [name for name in FREE_PARAMETERS if name not in normalised] + if missing: + raise KeyError( + f"The calibration parameters are missing {', '.join(missing)}; " + f"the decision vector needs {', '.join(FREE_PARAMETERS)}." + ) + return np.asarray([float(normalised[name]) for name in FREE_PARAMETERS], dtype=np.float64) + + +def is_admissible(vector: Sequence[float] | np.ndarray) -> bool: + """Whether a decision vector describes a configuration the model accepts. + + A vector is admissible when every free value lies inside its bound, when + the two searched weights add up to at most 1 and when the derived ``w_3`` + lies inside its own bound. The derived weight is taken from + :func:`vector_to_parameters`, the very value the run would be given, so + that the guard can never disagree with the derivation and let through a + candidate + :class:`rubem.configuration.calibration_parameters.CalibrationParameters` + would then refuse. + + :param vector: The eight free values, in the order of :data:`FREE_PARAMETERS`. + :type vector: collections.abc.Sequence[float] | numpy.ndarray + + :return: ``True`` when the model can be run with this candidate. + :rtype: bool + """ + values = np.asarray(vector, dtype=np.float64).ravel() + if values.size != len(FREE_PARAMETERS): + return False + if not np.all(np.isfinite(values)): + return False + for value, (minimum, maximum) in zip(values, bounds(), strict=True): + if not minimum <= value <= maximum: + return False + + weights = float(values[FREE_PARAMETERS.index("w_1")] + values[FREE_PARAMETERS.index("w_2")]) + if weights > 1.0: + return False + derived = vector_to_parameters(values)[DERIVED_PARAMETER] + minimum, maximum = variable_range(DERIVED_PARAMETER) + return minimum <= derived <= maximum + + +def weights_constraint() -> "LinearConstraint": + """Return the linear constraint ``w_1 + w_2 <= 1`` on the decision vector. + + The differential evolution receives it so that it never spends a model run + on a candidate whose slope factor weight would come out negative: SciPy + evaluates the constraint before the objective and skips the members that + violate it. + + :return: The constraint, over the eight entries of the decision vector. + :rtype: scipy.optimize.LinearConstraint + + :raises ImportError: If SciPy is not installed. + """ + from scipy.optimize import LinearConstraint + + coefficients = np.zeros((1, len(FREE_PARAMETERS)), dtype=np.float64) + coefficients[0, FREE_PARAMETERS.index("w_1")] = 1.0 + coefficients[0, FREE_PARAMETERS.index("w_2")] = 1.0 + return LinearConstraint(coefficients, -np.inf, 1.0) diff --git a/rubem/calibration/runner.py b/rubem/calibration/runner.py new file mode 100644 index 00000000..3d73fd21 --- /dev/null +++ b/rubem/calibration/runner.py @@ -0,0 +1,817 @@ +"""The calibration itself: a differential evolution over the free parameters. + +:func:`calibrate` loads and validates the configuration once, reads the observed +series, and hands SciPy's ``differential_evolution`` a process pool whose workers +run one model evaluation each and exit. The pool is the isolation layer: PCRaster +keeps its clone and its raster memory process-wide, so a worker that is reused +would accumulate the state of every run it made; ``spawn`` with +``max_tasks_per_child=1`` gives every evaluation a fresh interpreter. + +The run directory receives three artifacts: ``evaluations.csv``, one row per +evaluation, ``result.json``, the best candidate and the settings that found it, +and ``-calibrated.json``, the calibrated configuration in the +format the input was written in. The rows of the table are ordered by the +candidate they evaluated and not by the moment they were written, so that the +order of the rows does not depend on the order in which the workers happened to +finish. It does not make the table itself reproducible: the ``id``, ``pid`` and +``elapsed_seconds`` columns differ between two runs of one seed, and from the +first generation on the selection reads the objective values, so two runs +evaluate the same candidates only when the simulation itself is reproducible +(on the real basins that requires a fixed ``RASTERS.ldd``, since ``lddcreate`` +does not always derive the same directions twice; see the note on the LDD +raster in the user guide). + +Number of evaluations + The population of a ``sobol`` initialization is ``max(5, popsize * 8)`` + rounded up to the next power of two, since a Sobol' sequence is balanced + only over a power-of-two sample: with the default ``popsize`` of 15 that is + 128 members, so a calibration costs one model run per member and per + generation, ``128 * (maxiter + 1)`` runs at most. The number is logged when + the calibration starts so that the budget can be checked before it runs. + Members whose weights violate ``w_1 + w_2 <= 1`` are skipped by SciPy + without a model run, so the real count is usually lower. + +SciPy is an optional dependency, imported inside :func:`calibrate`; importing +this module without it works and raises nothing. +""" + +import csv +import functools +import json +import logging +import math +import multiprocessing +import os +import sys +import tempfile +from concurrent.futures import ProcessPoolExecutor +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any + +import numpy as np + +from .._paths import PathInput, as_path +from ..configuration.model_configuration_file_v1 import ( + VARIABLE_IDS, + Aggregation, + ModelConfigurationFileV1, +) +from ._worker import _INADMISSIBLE_ERROR, EvaluationContext, evaluate +from .objective import Series, read_series +from .parameters import ( + CALIBRATION_PARAMETERS, + FREE_PARAMETERS, + bounds, + parameters_to_vector, + vector_to_parameters, + weights_constraint, +) + +logger = logging.getLogger(__name__) + +EVALUATIONS_DIRNAME = "evaluations" +"""Name of the directory under the run directory that collects the JSON records.""" + +EVALUATIONS_CSV = "evaluations.csv" +RESULT_JSON = "result.json" + +EVALUATION_COLUMNS = ( + "id", + "pid", + *CALIBRATION_PARAMETERS, + "nse", + "objective", + "elapsed_seconds", + "error", +) +"""The columns of ``evaluations.csv``, one row per evaluation.""" + + +class CalibrationError(RuntimeError): + """A calibration that could not start, or that produced no usable evaluation.""" + + +@dataclass(frozen=True) +class CalibrationSettings: + """How the differential evolution is run. + + The differential evolution defaults follow the method the calibration of the + model was designed around; the remaining fields are the knobs the command + line exposes. + + :param variable: Id of the output variable compared with the observations. + ``arn``, the accumulated total runoff in m3/s, by default; ``rnf`` is + the total runoff per cell in mm. + :type variable: str + + :param spinup_steps: Number of initial time steps excluded from the + efficiency. + :type spinup_steps: int + + :param maxiter: Maximum number of generations. + :type maxiter: int + + :param popsize: Population multiplier; see the module docstring for the + number of members it produces. + :type popsize: int + + :param seed: Seed of the random generator of the search, for a reproducible + calibration. ``None`` leaves the generator unseeded. + :type seed: int | None + + :param workers: Number of worker processes. ``None``, the default, resolves + at call time to one less than the number of CPUs, at least one, which + leaves a core for the parent. + :type workers: int | None + + :param temp_dir: Directory the per-evaluation output directories are + created in. ``None`` uses the temporary directory of the system. + :type temp_dir: str | None + + :param init: Initial population: the name of a SciPy initialization + (``sobol`` by default) or an explicit ``(S, 8)`` array with ``S > 4``. + :type init: str | numpy.ndarray + + :param strategy: Differential evolution strategy. + :type strategy: str + + :param mutation: Mutation constant, a ``(min, max)`` pair for dithering. + :type mutation: tuple[float, float] + + :param recombination: Crossover probability. + :type recombination: float + + :param polish: Whether a local search refines the best candidate. ``False``: + the local search would spend model runs on gradients the objective, a + simulation read from a table, does not provide reliably. + :type polish: bool + """ + + variable: str = "arn" + spinup_steps: int = 0 + maxiter: int = 100 + popsize: int = 15 + seed: int | None = None + workers: int | None = None + temp_dir: str | None = None + init: str | np.ndarray = "sobol" + strategy: str = "best1exp" + mutation: tuple[float, float] = (0.5, 1.0) + recombination: float = 0.7 + polish: bool = False + + +@dataclass(frozen=True) +class CalibrationResult: + """What a finished calibration found and where it wrote it. + + :param best_parameters: The nine calibration parameters of the best + candidate, the derived ``w_3`` included. + :type best_parameters: dict[str, float] + + :param best_nse: The mean Nash-Sutcliffe efficiency of the best candidate, + ``None`` when no recorded evaluation reached one. + :type best_nse: float | None + + :param best_objective: The objective value of the best candidate. + :type best_objective: float + + :param evaluations: Number of objective evaluations the search spent. + :type evaluations: int + + :param generations: Number of generations the search ran. + :type generations: int + + :param success: Whether SciPy reports the search as successful. + :type success: bool + + :param message: The termination message of SciPy. + :type message: str + + :param run_dir: The run directory the artifacts were written to. + :type run_dir: pathlib.Path + + :param evaluations_csv: The table of every evaluation. + :type evaluations_csv: pathlib.Path + + :param result_json: The summary of the calibration. + :type result_json: pathlib.Path + + :param calibrated_config: The input configuration with the best parameters. + :type calibrated_config: pathlib.Path + """ + + best_parameters: dict[str, float] + best_nse: float | None + best_objective: float + evaluations: int + generations: int + success: bool + message: str + run_dir: Path + evaluations_csv: Path + result_json: Path + calibrated_config: Path + + +def _prepare_worker() -> None: + """Silence a worker process before it evaluates its candidate. + + The pool starts every worker as a fresh interpreter with no logging + configuration, so the warnings the loader emits for a configuration whose + inputs are not validated again would reach the terminal through the last + resort handler of :mod:`logging`, once per evaluation; the parent already + reported them when it validated the inputs. The PCRaster framework also + prints a newline when an interpreter exits, one per evaluation. A null + handler on the root logger and a standard output pointed at the null device + keep both out of the terminal: every failure of an evaluation is recorded in + its JSON record and summarized by the parent instead. + """ + logging.getLogger().addHandler(logging.NullHandler()) + sys.stdout = Path(os.devnull).open("w", encoding="utf-8") + + +def _report_failures(records: list[dict[str, Any]]) -> None: + """Log one summary of the evaluations that failed, if any did.""" + failures = [ + record["error"] + for record in records + if record.get("error") not in (None, _INADMISSIBLE_ERROR) + ] + if failures: + logger.warning( + "%d of %d evaluation(s) failed and were ranked behind every evaluated candidate; " + "the first failure was: %s. Every failure is recorded in the evaluations table.", + len(failures), + len(records), + failures[0], + ) + + +@dataclass +class _Progress: + """The per-generation callback: it only logs, it never stops the search.""" + + evaluations_dir: Path + generations: int = field(default=0) + + def __call__(self, intermediate_result) -> bool: + """Log the generation, the best objective so far and the records written. + + SciPy inspects the signature of the callback and passes the + intermediate result by keyword under exactly this name; returning + ``True`` would halt the search, so the callback always returns + ``False``. + """ + self.generations = int(getattr(intermediate_result, "nit", 0) or 0) + logger.info( + "Generation %d: best objective %.6g, %d evaluation(s) recorded.", + self.generations, + float(intermediate_result.fun), + len(_record_paths(self.evaluations_dir)), + ) + return False + + +def calibrate( + config_path: PathInput, + observed_path: PathInput, + run_dir: PathInput, + settings: CalibrationSettings = CalibrationSettings(), +) -> CalibrationResult: + """Calibrate a configuration against an observed series and write the artifacts. + + :param config_path: The configuration file to calibrate, legacy or format 1.0. + :type config_path: str | os.PathLike[str] + + :param observed_path: The observed series, in the CSV layout of the time + series tables of the model or as a PCRaster TSS file. + :type observed_path: str | os.PathLike[str] + + :param run_dir: Directory the artifacts of the calibration are written to. + It is created when it does not exist. + :type run_dir: str | os.PathLike[str] + + :param settings: How the search is run, defaults to :class:`CalibrationSettings`. + :type settings: CalibrationSettings, optional + + :return: The best candidate and the artifacts of the run. + :rtype: CalibrationResult + + :raises CalibrationError: If SciPy is not installed, if the configuration + cannot be calibrated as asked, if the observed series shares no time + step or no station with what the configuration will sample, or if every + evaluation of the run failed. + :raises ImportError: If PCRaster or GDAL are not installed. + :raises FileNotFoundError: If the configuration or the observed series is + not there. + :raises ConfigurationError: If the configuration carries blocking problems. + """ + differential_evolution = _require_scipy() + + from ..api import Model + + if settings.variable not in VARIABLE_IDS: + raise CalibrationError( + f"'{settings.variable}' is not an output variable of the model; " + f"calibrate against one of {', '.join(VARIABLE_IDS)}." + ) + if settings.spinup_steps < 0: + raise CalibrationError( + f"The number of spin-up steps must not be negative, got {settings.spinup_steps}." + ) + if settings.workers is not None and settings.workers < 1: + raise CalibrationError( + f"The calibration needs at least one worker process, got {settings.workers}." + ) + + config_file = as_path(config_path).absolute() + configuration = Model.from_file(config_file, validate_input=True).configuration + file_v1 = _as_v1(configuration) + _check_sample_locations(configuration, file_v1) + + observed = read_series(observed_path) + aggregation = file_v1.model_simulation_output.time_series_samples.aggregation + _check_observed_series( + observed, + first_step=configuration.simulation_period.first_step, + last_step=configuration.simulation_period.last_step, + spinup_steps=settings.spinup_steps, + aggregation=aggregation, + station_ids=( + None + if aggregation is Aggregation.ZONES + else _station_ids(configuration.raster_files.sample_locations) + ), + ) + x0 = parameters_to_vector(configuration.calibration_parameters.model_dump()) + + run_directory = as_path(run_dir).absolute() + evaluations_dir = run_directory / EVALUATIONS_DIRNAME + evaluations_dir.mkdir(parents=True, exist_ok=True) + previous = _record_paths(evaluations_dir) + if previous: + raise CalibrationError( + f"The run directory {run_directory} already holds {len(previous)} evaluation " + "record(s) of an earlier calibration; consolidating both runs into one table " + "would mix them. Choose an empty run directory." + ) + temp_dir = ( + as_path(settings.temp_dir).absolute() + if settings.temp_dir is not None + else Path(tempfile.gettempdir()) + ) + temp_dir.mkdir(parents=True, exist_ok=True) + + workers = ( + settings.workers if settings.workers is not None else max(1, (os.cpu_count() or 2) - 1) + ) + context = EvaluationContext( + document=file_v1.to_dict(), + base_dir=configuration.base_dir, + variable=settings.variable, + observed=observed, + spinup_steps=settings.spinup_steps, + temp_dir=str(temp_dir), + evaluations_dir=str(evaluations_dir), + validate_input=False, + ) + + planned_members = _population_size(settings) + logger.info( + "Calibrating %d free parameter(s) with %d population member(s) per generation and " + "at most %d generation(s): up to %d model run(s), on %d worker process(es).", + len(FREE_PARAMETERS), + planned_members, + settings.maxiter, + planned_members * (settings.maxiter + 1), + workers, + ) + + progress = _Progress(evaluations_dir) + executor = ProcessPoolExecutor( + max_workers=workers, + mp_context=multiprocessing.get_context("spawn"), + max_tasks_per_child=1, + initializer=_prepare_worker, + ) + try: + optimum = differential_evolution( + functools.partial(evaluate, context=context), + bounds(), + strategy=settings.strategy, + maxiter=settings.maxiter, + popsize=settings.popsize, + mutation=settings.mutation, + recombination=settings.recombination, + rng=settings.seed, + polish=settings.polish, + init=settings.init, + x0=x0, + updating="deferred", + workers=executor.map, + constraints=weights_constraint(), + callback=progress, + ) + finally: + executor.shutdown(wait=True) + + records = _read_records(evaluations_dir) + evaluations_csv = run_directory / EVALUATIONS_CSV + _write_evaluations_csv(evaluations_csv, records) + + _report_failures(records) + best = _best_record(records) + if best is None: + raise CalibrationError( + "Every evaluation of the calibration failed; the first failure was: " + f"{_first_error(records)}. The table of the evaluations is {evaluations_csv}." + ) + + # The size of the population SciPy actually built, not the size this module + # predicted for the budget: a named initialization resizes the population + # (``sobol`` rounds it up to a power of two) and the summary must describe + # the search that ran. + population_size = int(np.shape(optimum.population)[0]) + if population_size != planned_members: + logger.info( + "The search used %d population member(s), %d were announced.", + population_size, + planned_members, + ) + + best_parameters = vector_to_parameters(optimum.x) + best_objective = float(optimum.fun) + result_json = run_directory / RESULT_JSON + calibrated_config = run_directory / f"{config_file.stem}-calibrated.json" + _write_result_json( + result_json, + best_parameters=best_parameters, + best_nse=best["nse"], + best_objective=best_objective, + optimum=optimum, + settings=settings, + workers=workers, + temp_dir=temp_dir, + population_size=population_size, + ) + _write_calibrated_config(calibrated_config, configuration, file_v1, best_parameters) + + logger.info( + "Calibration finished after %d evaluation(s) in %d generation(s): best objective %.6g.", + int(optimum.nfev), + int(optimum.nit), + best_objective, + ) + return CalibrationResult( + best_parameters=best_parameters, + best_nse=best["nse"], + best_objective=best_objective, + evaluations=int(optimum.nfev), + generations=int(optimum.nit), + success=bool(optimum.success), + message=str(optimum.message), + run_dir=run_directory, + evaluations_csv=evaluations_csv, + result_json=result_json, + calibrated_config=calibrated_config, + ) + + +def _require_scipy(): + """Return ``scipy.optimize.differential_evolution``, or explain how to install it. + + :raises CalibrationError: If SciPy is not installed. + """ + from .. import _deps + + missing = _deps.missing_calibration_deps() + if missing: + raise CalibrationError(_deps.calibration_deps_message(missing)) + from scipy.optimize import differential_evolution + + return differential_evolution + + +def _as_v1(configuration) -> ModelConfigurationFileV1: + """Return the configuration as a format 1.0 model, converting a legacy one.""" + if configuration.file_v1 is not None: + return configuration.file_v1 + return ModelConfigurationFileV1.from_legacy(configuration.file) + + +def _check_sample_locations(configuration, file_v1: ModelConfigurationFileV1) -> None: + """Refuse a configuration whose time series would have nowhere to be sampled. + + :raises CalibrationError: If the raster the configured aggregation needs is + not part of the configuration. + """ + aggregation = file_v1.model_simulation_output.time_series_samples.aggregation + if aggregation is Aggregation.ZONES: + if not configuration.raster_files.zones: + raise CalibrationError( + "The configuration aggregates its time series over zones but names no " + "zones raster; the calibration has no station to compare." + ) + return + if not configuration.raster_files.sample_locations: + raise CalibrationError( + "The configuration names no sample locations raster; the calibration has no " + "station to compare the observed series with." + ) + + +def _station_ids(sample_locations: PathInput) -> set[str]: + """Return the station ids of a sample locations raster. + + The raster is a nominal map, so its cells read as integers. Neither the + no-data cells nor the zeroes are stations, which is what the validation of + the inputs already asks of a sample raster: its ids are ``1..N`` and its + background is missing or zero. Every remaining value becomes the string the + model writes in the header of its time series table, ``str(int(value))``. + + :param sample_locations: The sample locations raster of the configuration. + :type sample_locations: str | os.PathLike[str] + + :return: The distinct station ids, as strings. + :rtype: set[str] + """ + from osgeo import gdal + + gdal.UseExceptions() + gdal.AllRegister() + dataset = gdal.Open(str(as_path(sample_locations))) + try: + band = dataset.GetRasterBand(1) + values = np.unique(np.asarray(band.ReadAsArray())) + no_data = band.GetNoDataValue() + finally: + dataset = None + + if no_data is not None: + values = values[values != no_data] + return {str(int(value)) for value in values.tolist() if math.isfinite(value) and value != 0} + + +def _check_observed_series( + observed: Series, + *, + first_step: int, + last_step: int, + spinup_steps: int, + aggregation: Aggregation, + station_ids: set[str] | None, +) -> None: + """Refuse an observed series that cannot be compared with what the run samples. + + Every evaluation aligns the observed series with the table its own + simulation wrote, so a series that shares no step or no station with that + table makes every evaluation fail. The mismatch is the same for every + candidate, so it is settled here, in the parent, before a single model run + is spent on it. + + The compared steps are the simulated steps, ``first_step`` to ``last_step``, + that survive the spin-up. The stations are the ids the ``point`` and the + ``subcatchment`` aggregations sample, the values of the sample locations + raster; the ``zones`` aggregation remaps the zone ids to ``1..N`` at run + time (the mapping is written to ``zones_mapping.csv``), so the parent cannot + know the ids of its stations and the station check is skipped for it. + + :param observed: The observed series, as the parent read it. + :type observed: rubem.calibration.objective.Series + + :param first_step: Number of the first simulated step. + :type first_step: int + + :param last_step: Number of the last simulated step. + :type last_step: int + + :param spinup_steps: Number of initial time steps excluded from the + efficiency. + :type spinup_steps: int + + :param aggregation: The aggregation the time series of the run use. + :type aggregation: Aggregation + + :param station_ids: The ids of the stations the run will sample, or ``None`` + when they are not known before the run. + :type station_ids: set[str] | None + + :raises CalibrationError: If no observed step is compared, or if no observed + station id is one of the stations of the run. + """ + compared_steps = {step for step in range(first_step, last_step + 1) if step > spinup_steps} + observed_steps = {int(step) for step in np.asarray(observed.steps).tolist()} + if not compared_steps & observed_steps: + raise CalibrationError( + f"The observed series covers the step(s) {min(observed_steps)} to " + f"{max(observed_steps)}, and the configuration simulates the step(s) " + f"{first_step} to {last_step}, of which the {spinup_steps} spin-up step(s) " + "are excluded from the efficiency: the two series would have no time step " + "in common and every evaluation would fail." + ) + + if aggregation is Aggregation.ZONES or station_ids is None: + return + + observed_ids = set(observed.stations) + if not observed_ids & station_ids: + raise CalibrationError( + f"The observed series has the station(s) {', '.join(sorted(observed_ids)) or 'none'} " + f"and the configuration samples the station(s) " + f"{', '.join(sorted(station_ids)) or 'none'}: the two series would have no " + "station in common and every evaluation would fail." + ) + + foreign = sorted(observed_ids - station_ids) + if foreign: + logger.warning( + "The observed series has %d station(s) the configuration does not sample " + "(%s); they are ignored by the calibration.", + len(foreign), + ", ".join(foreign), + ) + + +def _population_size(settings: CalibrationSettings) -> int: + """Return the number of population members the search will use. + + The arithmetic mirrors SciPy's: an explicit population is used as it is, a + named initialization gives ``max(5, popsize * 8)`` members, and ``sobol`` + rounds that up to the next power of two. + """ + if not isinstance(settings.init, str): + return int(np.shape(settings.init)[0]) + members = max(5, settings.popsize * len(FREE_PARAMETERS)) + if settings.init == "sobol": + return 1 << (members - 1).bit_length() + return members + + +def _record_paths(evaluations_dir: Path) -> list[Path]: + """Return the JSON records written so far, without the files still being staged.""" + return sorted(path for path in evaluations_dir.glob("*.json") if not path.name.startswith(".")) + + +def _read_records(evaluations_dir: Path) -> list[dict[str, Any]]: + """Read every JSON record of the run, in the order of :func:`_record_order`. + + A record that cannot be read is reported and skipped: an unreadable record + of one evaluation must not hide the result of all the others. + """ + records = [] + for path in _record_paths(evaluations_dir): + try: + records.append(json.loads(path.read_text(encoding="utf-8"))) + except (OSError, ValueError) as error: + logger.warning("Skipping the unreadable evaluation record %s: %s", path, error) + records.sort(key=_record_order) + return records + + +def _record_order(record: dict[str, Any]) -> tuple[tuple[float, ...], str]: + """Return the sort key of a record: its parameters, then its identifier. + + The records are not ordered by the moment they were written: the workers run + in parallel, so that order changes from run to run. Ordering them by the + candidate they evaluated 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", ""), + *(_cell(parameters.get(name)) for name in CALIBRATION_PARAMETERS), + _cell(record.get("nse")), + _cell(record.get("objective")), + _cell(record.get("elapsed_seconds")), + record.get("error") or "", + ] + ) + + +def _cell(value) -> str: + """Return a numeric cell with full precision, or an empty cell for a missing value.""" + return "" if value is None else repr(float(value)) + + +def _best_record(records: list[dict[str, Any]]) -> dict[str, Any] | None: + """Return the record of the lowest objective among those that were evaluated.""" + evaluated = [record for record in records if not record.get("error")] + if not evaluated: + return None + return min(evaluated, key=lambda record: float(record["objective"])) + + +def _first_error(records: list[dict[str, Any]]) -> str: + """Return the error of the first record that carries one.""" + return next((str(record["error"]) for record in records if record.get("error")), "unknown") + + +def _write_result_json( + path: Path, + *, + best_parameters: dict[str, float], + best_nse: float | None, + best_objective: float, + optimum, + settings: CalibrationSettings, + workers: int, + temp_dir: Path, + population_size: int, +) -> None: + """Write the summary of the calibration.""" + document = { + "best_parameters": best_parameters, + "best_nse": best_nse, + "best_objective": best_objective, + "nfev": int(optimum.nfev), + "nit": int(optimum.nit), + "success": bool(optimum.success), + "message": str(optimum.message), + "population_size": population_size, + "settings": _settings_document(settings, workers, temp_dir), + } + path.write_text(json.dumps(document, indent=2) + "\n", encoding="utf-8") + + +def _settings_document( + settings: CalibrationSettings, workers: int, temp_dir: Path +) -> dict[str, Any]: + """Return the settings as JSON-able data, resolved the way the run used them. + + The number of workers and the temporary directory are the values the run + actually used, not the ``None`` the caller may have left them at, so that + the summary describes a calibration that can be repeated. + """ + init = settings.init if isinstance(settings.init, str) else np.asarray(settings.init).tolist() + return { + "variable": settings.variable, + "spinup_steps": settings.spinup_steps, + "maxiter": settings.maxiter, + "popsize": settings.popsize, + "seed": settings.seed, + "workers": workers, + "temp_dir": str(temp_dir), + "init": init, + "strategy": settings.strategy, + "mutation": list(settings.mutation), + "recombination": settings.recombination, + "polish": settings.polish, + } + + +def _write_calibrated_config( + path: Path, + configuration, + file_v1: ModelConfigurationFileV1, + best_parameters: dict[str, float], +) -> None: + """Write the calibrated configuration in the format the input was written in. + + The document goes through the configuration models, so what is written is a + configuration the loader accepts. Its paths are the ones the loader resolved, + which are absolute even when the input file wrote them relative to its own + directory. + """ + calibrated = file_v1.model_copy( + update={ + "model_calibration_parameters": type(file_v1.model_calibration_parameters)( + alpha=best_parameters["alpha"], + b=best_parameters["beta"], + w_1=best_parameters["w_1"], + w_2=best_parameters["w_2"], + w_3=best_parameters["w_3"], + rcd=best_parameters["rcd"], + f=best_parameters["f"], + alpha_gw=best_parameters["alpha_gw"], + x=best_parameters["x"], + ) + } + ) + document = ( + calibrated.to_dict() + if configuration.file_v1 is not None + else calibrated.to_legacy().to_dict() + ) + path.write_text(json.dumps(document, indent=2) + "\n", encoding="utf-8") diff --git a/rubem/cli.py b/rubem/cli.py index fb9eb5fe..a1e0044e 100644 --- a/rubem/cli.py +++ b/rubem/cli.py @@ -1,6 +1,7 @@ """The ``rubem`` command line. -``rubem run -c [-s]`` runs a simulation; ``rubem config schema`` +``rubem run -c [-s]`` runs a simulation; ``rubem calibrate`` +fits the calibration parameters to an observed series; ``rubem config schema`` prints the JSON Schema of the configuration file. The former ``rubem -c `` spelling still works for one minor release and emits a ``DeprecationWarning``. @@ -21,7 +22,7 @@ from pydantic import ValidationError from . import __release__ -from ._deps import require_runtime_deps +from ._deps import require_calibration_deps, require_runtime_deps from .configuration._problems import ConfigurationError from .configuration.app_settings import AppSettings from .preprocessing.cli import app as preprocess_app @@ -178,6 +179,113 @@ def run( logger.info("RUBEM successfully finished!") +@app.command() +def calibrate( + configfile: Annotated[ + Path, + typer.Option( + "-c", + "--configfile", + callback=_configfile_callback, + help="Path to the configuration file (JSON).", + ), + ], + observed: Annotated[ + Path, + typer.Option( + "--observed", + exists=True, + dir_okay=False, + help="Observed series at the sample stations (CSV or PCRaster TSS).", + ), + ], + run_dir: Annotated[ + Path, + typer.Option("-o", "--run-dir", help="Directory for the calibration artifacts."), + ], + variable: Annotated[ + str, typer.Option("--variable", help="Output variable compared with the observed series.") + ] = "arn", + spinup_steps: Annotated[ + int, + typer.Option("--spinup-steps", min=0, help="Leading time steps excluded from the NSE."), + ] = 0, + maxiter: Annotated[ + int, typer.Option("--maxiter", min=1, help="Maximum number of generations.") + ] = 100, + popsize: Annotated[ + int, typer.Option("--popsize", min=1, help="Population multiplier of the search.") + ] = 15, + seed: Annotated[ + int | None, typer.Option("--seed", help="Seed of the differential evolution.") + ] = None, + workers: Annotated[ + int | None, + typer.Option("--workers", min=1, help="Parallel evaluations (default: all cores but one)."), + ] = None, + temp_dir: Annotated[ + Path | None, + typer.Option("--temp-dir", help="Parent of the per-evaluation output directories."), + ] = None, +) -> None: + """Calibrate the model parameters against an observed series.""" + require_runtime_deps() + require_calibration_deps() + + from .calibration.runner import CalibrationError, CalibrationSettings + from .calibration.runner import calibrate as run_calibration + + settings_arguments: dict = { + "variable": variable, + "spinup_steps": spinup_steps, + "maxiter": maxiter, + "popsize": popsize, + "seed": seed, + } + if workers is not None: + settings_arguments["workers"] = workers + if temp_dir is not None: + settings_arguments["temp_dir"] = str(temp_dir) + + try: + # The library only logs; the progress a command-line calibration is + # expected to show is written here, so that an embedded calibration + # stays silent unless its host configures logging. + print("Loading configuration and validating inputs...", flush=True) + print("Calibration started...", flush=True) + result = run_calibration( + configfile, observed, run_dir, CalibrationSettings(**settings_arguments) + ) + except (ConfigurationError, ValidationError, ValueError) as e: + # A configuration the user can fix: no traceback, the message says what. + logger.critical("Invalid configuration: %s", e) + raise typer.Exit(code=1) from e + except CalibrationError as e: + # The calibration itself refused to run or produced nothing usable; + # the message says what, so no traceback either. + logger.critical("Calibration failed: %s", e) + raise typer.Exit(code=1) from e + except KeyboardInterrupt as e: + logger.critical("RUBEM was interrupted by the user.") + raise typer.Exit(code=2) from e + except Exception as e: + logger.critical("RUBEM unexpectedly quit.") + logger.exception(e) + raise typer.Exit(code=1) from e + + best_nse = "n/a" if result.best_nse is None else f"{result.best_nse:.6f}" + print("Calibration finished successfully!", flush=True) + print(f"Best NSE: {best_nse}", flush=True) + print(f"Best objective: {result.best_objective:.6g}", flush=True) + print(f"Evaluations: {result.evaluations}", flush=True) + print(f"Generations: {result.generations}", flush=True) + print(f"Evaluations table: {result.evaluations_csv}", flush=True) + print(f"Result: {result.result_json}", flush=True) + print(f"Calibrated configuration: {result.calibrated_config}", flush=True) + + logger.info("RUBEM successfully finished!") + + @config_app.command("schema") def config_schema( schema_format: Annotated[ @@ -258,7 +366,7 @@ def _legacy_arguments(argv: Sequence[str]) -> list[str]: global _legacy_invocation _legacy_invocation = False arguments = list(argv) - if not arguments or arguments[0] in ("run", "config"): + if not arguments or arguments[0] in ("run", "config", "calibrate", "preprocess"): return arguments options = [a for a in arguments if a.startswith("-")] if not options or not {o.split("=", 1)[0] for o in options} <= _LEGACY_OPTIONS: diff --git a/tests/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..bc76f666 --- /dev/null +++ b/tests/unit/calibration/test_objective.py @@ -0,0 +1,396 @@ +"""The calibration objective: masks, Nash-Sutcliffe efficiency, tables and alignment.""" + +import numpy as np +import pytest + +from rubem.calibration.objective import ( + INADMISSIBLE_OBJECTIVE, + MISSING_VALUES, + PCRASTER_MISSING, + Series, + evaluate_series, + nash_sutcliffe, + objective, + read_series, + valid_mask, +) + +# The PCRaster missing value as it comes back from a Float32 raster or table: +# 1e31 has no exact Float32 representation, so what is read is slightly below it. +FLOAT32_MISSING = float(np.float32(1e31)) + +# Every gap marker, used wherever a test drops a pair through one of them. +GAPS = (float("nan"), -9999.0, 1e31, FLOAT32_MISSING) + + +def series(steps, **stations): + """A :class:`Series` built from plain lists, as the readers produce it.""" + return Series( + steps=np.asarray(steps, dtype=np.int64), + stations={ + station: np.asarray(values, dtype=np.float64) for station, values in stations.items() + }, + ) + + +class TestValidMask: + @pytest.mark.unit + def test_ordinary_values_are_valid(self): + assert valid_mask(np.array([0.0, -1.5, 3.25, 1e29])).tolist() == [True] * 4 + + @pytest.mark.unit + @pytest.mark.parametrize("gap", GAPS) + def test_every_gap_marker_is_masked(self, gap): + assert valid_mask(np.array([1.0, gap, 2.0])).tolist() == [True, False, True] + + @pytest.mark.unit + def test_the_documented_markers_are_the_ones_the_mask_rejects(self): + assert MISSING_VALUES == (-9999.0,) + assert PCRASTER_MISSING == 1e31 + assert not valid_mask(np.array([PCRASTER_MISSING])).any() + + @pytest.mark.unit + def test_infinities_are_masked(self): + assert valid_mask(np.array([np.inf, -np.inf])).tolist() == [False, False] + + @pytest.mark.unit + def test_the_threshold_itself_is_the_first_value_rejected(self): + # The mask keeps everything strictly below 1e30 and drops 1e30 itself, + # which is what makes the Float32 round-trip of 1e31 a gap as well. + below = float(np.nextafter(1e30, 0.0)) + + assert valid_mask(np.array([below, 1e30])).tolist() == [True, False] + assert below < 1e30 < FLOAT32_MISSING + + +class TestNashSutcliffe: + @pytest.mark.unit + def test_an_exact_simulation_gives_one(self): + # residuals 0 / variability 2 -> NSE = 1 - 0 = 1.0 + observed = np.array([1.0, 2.0, 3.0]) + + assert nash_sutcliffe(observed.copy(), observed) == 1.0 + + @pytest.mark.unit + def test_the_mean_of_the_observations_gives_zero(self): + # observed [1, 2, 3], mean 2; simulated [2, 2, 2] + # residuals (1 + 0 + 1) = 2, variability (1 + 0 + 1) = 2 -> NSE = 1 - 2/2 = 0.0 + assert nash_sutcliffe(np.array([2.0, 2.0, 2.0]), np.array([1.0, 2.0, 3.0])) == 0.0 + + @pytest.mark.unit + def test_a_reversed_simulation_gives_a_negative_efficiency(self): + # observed [1, 2, 3], mean 2; simulated [3, 2, 1] + # residuals (4 + 0 + 4) = 8, variability 2 -> NSE = 1 - 8/2 = -3.0 + assert nash_sutcliffe(np.array([3.0, 2.0, 1.0]), np.array([1.0, 2.0, 3.0])) == -3.0 + + @pytest.mark.unit + @pytest.mark.parametrize("gap", GAPS) + def test_a_gap_in_the_observations_drops_the_pair(self, gap): + # The fourth pair is dropped, the first three are exact -> NSE = 1.0 + simulated = np.array([1.0, 2.0, 3.0, 99.0]) + observed = np.array([1.0, 2.0, 3.0, gap]) + + assert nash_sutcliffe(simulated, observed) == 1.0 + + @pytest.mark.unit + @pytest.mark.parametrize("gap", GAPS) + def test_a_gap_in_the_simulation_drops_the_pair(self, gap): + # Dropping the fourth pair leaves observed [1, 2, 3] with mean 2 and + # simulated [3, 2, 1]: residuals 8, variability 2 -> NSE = -3.0 + simulated = np.array([3.0, 2.0, 1.0, gap]) + observed = np.array([1.0, 2.0, 3.0, 4.0]) + + assert nash_sutcliffe(simulated, observed) == -3.0 + + @pytest.mark.unit + def test_a_gap_is_dropped_from_the_mean_of_the_observations(self): + # The third pair is dropped by the simulation, so the mean is that of + # [1, 3] = 2 and not that of [1, 3, 8]: variability (1 + 1) = 2, + # residuals (1 + 1) = 2 -> NSE = 1 - 2/2 = 0.0 + simulated = np.array([2.0, 2.0, float("nan")]) + observed = np.array([1.0, 3.0, 8.0]) + + assert nash_sutcliffe(simulated, observed) == 0.0 + + @pytest.mark.unit + def test_a_single_valid_pair_has_no_efficiency(self): + assert nash_sutcliffe(np.array([1.0, 2.0]), np.array([1.0, float("nan")])) is None + + @pytest.mark.unit + def test_a_station_without_one_valid_pair_has_no_efficiency(self): + # No pair survives the mask at all: the count guard must answer before + # anything is read out of the empty selection. + assert nash_sutcliffe(np.array([1.0, -9999.0]), np.array([np.nan, 2.0])) is None + assert nash_sutcliffe(np.array([np.nan, np.nan]), np.array([np.nan, np.nan])) is None + + @pytest.mark.unit + def test_two_valid_pairs_are_enough(self): + # observed [1, 3], mean 2, variability 2; simulated [1, 3] -> NSE = 1.0 + simulated = np.array([1.0, 3.0, 5.0]) + observed = np.array([1.0, 3.0, -9999.0]) + + assert nash_sutcliffe(simulated, observed) == 1.0 + + @pytest.mark.unit + def test_observations_without_variance_have_no_efficiency(self): + # The sum of the squared deviations of [0.1, 0.1, 0.1] is not exactly + # zero in floating point, so constancy is what the guard must test. + assert nash_sutcliffe(np.array([0.1, 0.2, 0.3]), np.array([0.1, 0.1, 0.1])) is None + + @pytest.mark.unit + def test_only_the_valid_observations_decide_the_variance(self): + # The valid observations are [5, 5]: constant, so there is no efficiency + # although the series as written also carries a 7. + assert nash_sutcliffe(np.array([5.0, 5.0, np.nan]), np.array([5.0, 5.0, 7.0])) is None + + @pytest.mark.unit + def test_series_of_different_lengths_are_refused(self): + with pytest.raises(ValueError, match="same shape"): + nash_sutcliffe(np.array([1.0, 2.0]), np.array([1.0, 2.0, 3.0])) + + +class TestObjective: + @pytest.mark.unit + @pytest.mark.parametrize( + ("nse", "expected"), + [ + # 1000 * (100 * (1 - 1))^2 = 0 + (1.0, 0.0), + # 1000 * (100 * 1)^2 = 1000 * 10000 = 1e7 + (0.0, 1.0e7), + # 1000 * (100 * 2)^2 = 1000 * 40000 = 4e7 + (-1.0, 4.0e7), + ], + ) + def test_the_objective_of_a_known_efficiency(self, nse, expected): + assert objective(nse) == expected + + @pytest.mark.unit + def test_the_inadmissible_objective_is_far_above_any_run(self): + # An NSE of -100 is an absurd simulation and still gives ~1.02e11. + assert objective(-100.0) < INADMISSIBLE_OBJECTIVE / 1e10 + assert np.isfinite(INADMISSIBLE_OBJECTIVE) + + +class TestReadSeries: + @pytest.mark.unit + def test_a_handwritten_table_in_the_model_csv_layout(self, tmp_path): + path = tmp_path / "tss_arn.csv" + path.write_text("0;1;2\n1;1.5;2.5\n2;3.5;4.5\n3;5.5;6.5\n", encoding="utf8") + + table = read_series(path) + + assert table.steps.tolist() == [1, 2, 3] + assert sorted(table.stations) == ["1", "2"] + assert table.stations["1"].tolist() == [1.5, 3.5, 5.5] + assert table.stations["2"].tolist() == [2.5, 4.5, 6.5] + + @pytest.mark.unit + def test_gaps_of_a_csv_table_are_kept_as_written(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;A;B\n1;-9999;2.0\n2;1e31;\n", encoding="utf8") + + table = read_series(path) + + assert table.stations["A"].tolist() == [-9999.0, 1e31] + assert table.stations["B"][0] == 2.0 + assert np.isnan(table.stations["B"][1]), "an empty cell is a gap" + assert valid_mask(table.stations["A"]).tolist() == [False, False] + + @pytest.mark.unit + def test_the_ids_of_a_real_run_are_the_ids_of_the_sample_raster(self, tmp_path): + """The table the model writes must be readable exactly as it lands.""" + from tests.unit.core.test_core import run_model + + run_model(str(tmp_path)) + + table = read_series(tmp_path / "out" / "tss_arn.csv") + + # The synthetic dataset marks two sample locations, with the ids 1 and 2. + assert sorted(table.stations) == ["1", "2"] + assert table.steps.tolist() == [1, 2] + for values in table.stations.values(): + assert values.shape == (2,) + assert valid_mask(values).all() + + @pytest.mark.unit + def test_the_windows_line_endings_the_model_writes(self, tmp_path): + # rubem.file._file_conversions.tss2csv writes through csv.writer, whose + # line terminator is CRLF; the reader must not take the carriage return + # for part of the last value of a row. + path = tmp_path / "tss_arn.csv" + path.write_bytes(b"0;1;2\r\n1;1.5;2.5\r\n2;3.5;4.5\r\n") + + table = read_series(path) + + assert table.steps.tolist() == [1, 2] + assert table.stations["2"].tolist() == [2.5, 4.5] + + @pytest.mark.unit + def test_the_exponent_spellings_of_the_pcraster_missing_value(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;A\n1;1e+31\n2;1E31\n3;2.0\n", encoding="utf8") + + table = read_series(path) + + assert valid_mask(table.stations["A"]).tolist() == [False, False, True] + + @pytest.mark.unit + def test_a_handwritten_pcraster_time_series(self, tmp_path): + path = tmp_path / "observed.tss" + path.write_text( + "Observed streamflow\n3\ntimestep\n11\n22\n1 1.5 2.5\n2 3.5 4.5\n", + encoding="utf8", + ) + + table = read_series(path) + + assert table.steps.tolist() == [1, 2] + assert sorted(table.stations) == ["11", "22"] + assert table.stations["11"].tolist() == [1.5, 3.5] + assert table.stations["22"].tolist() == [2.5, 4.5] + + @pytest.mark.unit + def test_a_pcraster_time_series_without_the_header(self, tmp_path): + # The form the model itself writes: every line is a data row and the + # stations are numbered in column order. + path = tmp_path / "observed.tss" + path.write_text("1 1.5 2.5\n2 3.5 4.5\n", encoding="utf8") + + table = read_series(path) + + assert table.steps.tolist() == [1, 2] + assert sorted(table.stations) == ["1", "2"] + assert table.stations["2"].tolist() == [2.5, 4.5] + + @pytest.mark.unit + def test_a_row_that_does_not_start_with_a_time_step_names_the_file(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;A\n2000-01-01;1.0\n2000-01-02;2.0\n", encoding="utf8") + + with pytest.raises(ValueError, match="Row 2 .* does not start with a time step"): + read_series(path) + + @pytest.mark.unit + def test_a_repeated_station_id_is_refused(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;1;1\n1;1.0;2.0\n", encoding="utf8") + + with pytest.raises(ValueError, match="repeated station ids: 1"): + read_series(path) + + @pytest.mark.unit + def test_a_row_that_does_not_match_the_header_is_refused(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;1;2\n1;1.0;2.0\n2;3.0\n", encoding="utf8") + + with pytest.raises(ValueError, match="Row 3"): + read_series(path) + + @pytest.mark.unit + def test_an_empty_table_is_refused(self, tmp_path): + path = tmp_path / "observed.csv" + path.write_text("", encoding="utf8") + + with pytest.raises(ValueError, match="empty"): + read_series(path) + + +class TestEvaluateSeries: + @pytest.mark.unit + def test_the_mean_over_the_stations(self): + # Station A: exact -> 1.0. Station B: observed [1, 2, 3] with the + # simulation reversed -> -3.0. Mean = (1.0 - 3.0) / 2 = -1.0 + simulated = series([1, 2, 3], A=[4.0, 5.0, 6.0], B=[3.0, 2.0, 1.0]) + observed = series([1, 2, 3], A=[4.0, 5.0, 6.0], B=[1.0, 2.0, 3.0]) + + mean, per_station = evaluate_series(simulated, observed) + + assert per_station == {"A": 1.0, "B": -3.0} + assert mean == -1.0 + + @pytest.mark.unit + def test_only_the_shared_steps_are_compared(self): + # The shared steps are 3 and 4: observed [1, 3] with mean 2 and + # variability 2, simulated [1, 3] -> NSE = 1.0 + simulated = series([1, 2, 3, 4], A=[99.0, 99.0, 1.0, 3.0]) + observed = series([3, 4, 5], A=[1.0, 3.0, 99.0]) + + mean, per_station = evaluate_series(simulated, observed) + + assert per_station == {"A": 1.0} + assert mean == 1.0 + + @pytest.mark.unit + def test_the_spin_up_steps_are_excluded(self): + # Without the spin-up the first step would make the simulation exact; + # with two spin-up steps only the steps 3 and 4 are compared: + # observed [1, 3], mean 2, variability 2; simulated [3, 1] + # -> residuals (4 + 4) = 8 -> NSE = 1 - 8/2 = -3.0 + simulated = series([1, 2, 3, 4], A=[0.0, 0.0, 3.0, 1.0]) + observed = series([1, 2, 3, 4], A=[0.0, 0.0, 1.0, 3.0]) + + mean, per_station = evaluate_series(simulated, observed, spinup_steps=2) + + assert per_station == {"A": -3.0} + assert mean == -3.0 + + @pytest.mark.unit + def test_a_station_present_on_one_side_only_is_ignored(self): + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[1.0, 2.0, 3.0]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0], C=[1.0, 2.0, 3.0]) + + mean, per_station = evaluate_series(simulated, observed) + + assert per_station == {"A": 1.0} + assert mean == 1.0 + + @pytest.mark.unit + def test_a_station_without_variance_is_reported_and_left_out_of_the_mean(self): + # Station A is exact -> 1.0; station B has constant observations -> None. + # The mean is that of the stations with a value: 1.0 + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[0.1, 0.2, 0.3]) + observed = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[0.1, 0.1, 0.1]) + + mean, per_station = evaluate_series(simulated, observed) + + assert per_station == {"A": 1.0, "B": None} + assert mean == 1.0 + + @pytest.mark.unit + def test_no_shared_station_names_both_sides(self): + simulated = series([1, 2], A=[1.0, 2.0]) + observed = series([1, 2], B=[1.0, 2.0]) + + with pytest.raises(ValueError, match="no station in common") as error: + evaluate_series(simulated, observed) + + assert "A" in str(error.value) and "B" in str(error.value) + + @pytest.mark.unit + def test_no_shared_step_names_the_steps(self): + simulated = series([1, 2], A=[1.0, 2.0]) + observed = series([7, 8], A=[1.0, 2.0]) + + with pytest.raises(ValueError, match="no time step in common") as error: + evaluate_series(simulated, observed) + + assert "[1, 2]" in str(error.value) and "[7, 8]" in str(error.value) + + @pytest.mark.unit + def test_a_spin_up_that_eats_every_step_is_an_error(self): + simulated = series([1, 2], A=[1.0, 2.0]) + observed = series([1, 2], A=[1.0, 2.0]) + + with pytest.raises(ValueError, match="no time step in common"): + evaluate_series(simulated, observed, spinup_steps=2) + + @pytest.mark.unit + def test_no_station_with_an_efficiency_is_an_error(self): + simulated = series([1, 2, 3], A=[1.0, 2.0, 3.0], B=[1.0, 2.0, 3.0]) + observed = series([1, 2, 3], A=[5.0, 5.0, 5.0], B=[1.0, -9999.0, np.nan]) + + with pytest.raises(ValueError, match="No station has a Nash-Sutcliffe") as error: + evaluate_series(simulated, observed) + + assert "A, B" in str(error.value) diff --git a/tests/unit/calibration/test_parameters.py b/tests/unit/calibration/test_parameters.py new file mode 100644 index 00000000..90d0db87 --- /dev/null +++ b/tests/unit/calibration/test_parameters.py @@ -0,0 +1,227 @@ +"""The decision vector of the calibration: bounds, conversions and admissibility.""" + +import numpy as np +import pytest + +from rubem.calibration.parameters import ( + CALIBRATION_PARAMETERS, + DERIVED_PARAMETER, + FREE_PARAMETERS, + bounds, + is_admissible, + parameters_to_vector, + vector_to_parameters, + weights_constraint, +) +from rubem.configuration.app_settings import AppSettings +from rubem.configuration.calibration_parameters import CalibrationParameters +from rubem.configuration.model_configuration import ModelConfiguration +from tests.helpers.synthetic import write_synthetic_dataset + +# The ranges of rubem/appsettings.json, written out so that a change of the +# settings file, or of the order of the decision vector, has to be made here +# too instead of passing unnoticed through a comparison with the same source. +SETTINGS_BOUNDS = { + "alpha": (0.01, 10.0), + "beta": (0.01, 1.0), + "w_1": (0.0, 1.0), + "w_2": (0.0, 1.0), + "w_3": (0.0, 1.0), + "rcd": (1.0, 10.0), + "f": (0.01, 1.0), + "alpha_gw": (0.01, 1.0), + "x": (0.0, 1.0), +} + + +def admissible_vector(): + """A decision vector well inside every bound, with weights that add up below 1.""" + return np.array([4.5, 0.5, 0.333, 0.333, 5.0, 0.5, 0.5, 0.5], dtype=np.float64) + + +class TestBounds: + @pytest.mark.unit + def test_the_bounds_are_the_ranges_of_the_application_settings(self): + ranges = AppSettings.default().value_ranges.variables + + assert bounds() == [(ranges[name].min, ranges[name].max) for name in FREE_PARAMETERS] + + @pytest.mark.unit + def test_the_bounds_are_the_ones_written_in_the_settings_file(self): + assert bounds() == [SETTINGS_BOUNDS[name] for name in FREE_PARAMETERS] + assert [name for name in SETTINGS_BOUNDS if name != DERIVED_PARAMETER] == list( + FREE_PARAMETERS + ) + + @pytest.mark.unit + def test_the_free_parameters_do_not_include_the_derived_weight(self): + assert DERIVED_PARAMETER not in FREE_PARAMETERS + assert set(CALIBRATION_PARAMETERS) == set(FREE_PARAMETERS) | {DERIVED_PARAMETER} + assert len(bounds()) == len(FREE_PARAMETERS) == 8 + + +class TestConversions: + @pytest.mark.unit + def test_a_vector_becomes_the_nine_parameters_with_the_derived_weight(self): + parameters = vector_to_parameters(admissible_vector()) + + assert list(parameters) == list(CALIBRATION_PARAMETERS) + assert parameters["beta"] == 0.5 + assert parameters["w_3"] == 1.0 - (0.333 + 0.333) + assert parameters["w_1"] + parameters["w_2"] + parameters["w_3"] == pytest.approx(1.0) + + @pytest.mark.unit + def test_a_vector_of_the_wrong_size_is_rejected(self): + with pytest.raises(ValueError, match="must have 8 entries"): + vector_to_parameters([1.0, 2.0, 3.0]) + + @pytest.mark.unit + def test_the_configuration_parameters_round_trip_through_the_vector(self, tmp_path): + config = write_synthetic_dataset(str(tmp_path)) + configuration = ModelConfiguration(config, validate_input=False) + parameters = configuration.calibration_parameters.model_dump() + + vector = parameters_to_vector(parameters) + + assert vector_to_parameters(vector) == pytest.approx(parameters, abs=1e-12) + + @pytest.mark.unit + def test_the_legacy_spellings_of_the_parameters_are_accepted(self): + legacy = { + "alpha": 4.5, + "b": 0.5, + "w1": 0.333, + "w2": 0.333, + "w3": 0.334, + "rcd": 5.0, + "f": 0.5, + "alpha_gw": 0.5, + "x": 0.5, + } + + assert parameters_to_vector(legacy) == pytest.approx(admissible_vector()) + + @pytest.mark.unit + def test_a_missing_parameter_names_what_the_vector_needs(self): + with pytest.raises(KeyError, match="alpha_gw"): + parameters_to_vector({"alpha": 4.5, "beta": 0.5}) + + +class TestAdmissibility: + @pytest.mark.unit + def test_a_vector_inside_the_bounds_is_admissible(self): + assert is_admissible(admissible_vector()) + + @pytest.mark.unit + def test_a_vector_on_the_bounds_is_admissible(self): + vector = np.array([minimum for minimum, _ in bounds()], dtype=np.float64) + + assert is_admissible(vector) + + @pytest.mark.unit + @pytest.mark.parametrize("index", range(8)) + def test_a_value_above_its_bound_is_rejected(self, index): + vector = admissible_vector() + vector[index] = bounds()[index][1] + 1.0 + + assert not is_admissible(vector) + + @pytest.mark.unit + @pytest.mark.parametrize("index", range(8)) + def test_a_value_below_its_bound_is_rejected(self, index): + vector = admissible_vector() + vector[index] = bounds()[index][0] - 1.0 + + assert not is_admissible(vector) + + @pytest.mark.unit + def test_weights_that_add_up_above_one_are_rejected(self): + vector = admissible_vector() + vector[FREE_PARAMETERS.index("w_1")] = 0.6 + vector[FREE_PARAMETERS.index("w_2")] = 0.6 + + # Both weights are inside their own [0, 1] bound; what makes the + # candidate inadmissible is their sum, and the negative w_3 it derives. + assert vector_to_parameters(vector)["w_3"] < 0.0 + assert not is_admissible(vector) + + @pytest.mark.unit + def test_weights_that_add_up_to_exactly_one_are_admissible(self): + vector = admissible_vector() + vector[FREE_PARAMETERS.index("w_1")] = 0.4 + vector[FREE_PARAMETERS.index("w_2")] = 0.6 + + assert vector_to_parameters(vector)["w_3"] == 0.0 + assert is_admissible(vector) + + @pytest.mark.unit + @pytest.mark.parametrize("hundredths", range(101)) + def test_every_pair_of_weights_that_adds_up_to_one_reaches_the_model(self, hundredths): + """What the guard admits, the configuration model must accept. + + ``1 - w_1 - w_2`` rounds twice and lands one ulp below zero for 20 of + these 101 pairs, ``w_1 = 0.33, w_2 = 0.67`` among them, which is a + weight :class:`CalibrationParameters` refuses; deriving the weight with + a single subtraction gives exactly zero for all of them. + """ + vector = admissible_vector() + vector[FREE_PARAMETERS.index("w_1")] = hundredths / 100.0 + vector[FREE_PARAMETERS.index("w_2")] = (100 - hundredths) / 100.0 + parameters = vector_to_parameters(vector) + + assert parameters["w_3"] == 0.0 + assert is_admissible(vector) + # No exception: the candidate the guard admits is one the model runs. + CalibrationParameters(**parameters) + + @pytest.mark.unit + def test_the_guard_and_the_derivation_cannot_disagree(self): + """Every admissible candidate of a grid is one the model accepts.""" + refused = [] + for first in range(0, 101, 7): + for second in range(0, 101, 11): + vector = admissible_vector() + vector[FREE_PARAMETERS.index("w_1")] = first / 100.0 + vector[FREE_PARAMETERS.index("w_2")] = second / 100.0 + if not is_admissible(vector): + continue + try: + CalibrationParameters(**vector_to_parameters(vector)) + except ValueError: + refused.append((first / 100.0, second / 100.0)) + + assert not refused, f"admitted but refused by the model: {refused}" + + @pytest.mark.unit + def test_a_vector_of_the_wrong_size_or_with_a_nan_is_rejected(self): + vector = admissible_vector() + vector[0] = np.nan + + assert not is_admissible(vector) + assert not is_admissible(admissible_vector()[:-1]) + + +class TestWeightsConstraint: + """The only part of the module that needs the optional SciPy extra.""" + + @pytest.fixture(autouse=True) + def _scipy(self): + pytest.importorskip("scipy.optimize") + + @pytest.mark.unit + def test_the_constraint_is_the_sum_of_the_two_searched_weights(self): + constraint = weights_constraint() + + expected = np.zeros((1, len(FREE_PARAMETERS))) + expected[0, FREE_PARAMETERS.index("w_1")] = 1.0 + expected[0, FREE_PARAMETERS.index("w_2")] = 1.0 + assert np.array_equal(np.asarray(constraint.A), expected) + assert np.all(np.asarray(constraint.lb) == -np.inf) + assert np.all(np.asarray(constraint.ub) == 1.0) + + @pytest.mark.unit + def test_the_constraint_evaluates_to_the_sum_of_the_weights(self): + constraint = weights_constraint() + vector = admissible_vector() + + assert np.asarray(constraint.A) @ vector == pytest.approx([0.666]) diff --git a/tests/unit/calibration/test_runner.py b/tests/unit/calibration/test_runner.py new file mode 100644 index 00000000..0d340619 --- /dev/null +++ b/tests/unit/calibration/test_runner.py @@ -0,0 +1,729 @@ +"""The calibration runner: a small deterministic differential evolution and its artifacts.""" + +import csv +import inspect +import json +import logging +import multiprocessing +import os +import shutil +import subprocess +import sys +from concurrent.futures import ProcessPoolExecutor +from pathlib import Path + +import numpy as np +import pytest + +from rubem import _deps +from rubem.api import Model +from rubem.calibration import runner +from rubem.calibration.objective import Series +from rubem.calibration.parameters import ( + CALIBRATION_PARAMETERS, + FREE_PARAMETERS, + bounds, + parameters_to_vector, +) +from rubem.calibration.runner import ( + EVALUATION_COLUMNS, + CalibrationError, + CalibrationSettings, + _check_observed_series, + _population_size, + _Progress, + _station_ids, + calibrate, +) +from rubem.configuration.model_configuration import ModelConfiguration +from rubem.configuration.model_configuration_file import ModelConfigurationFile +from rubem.configuration.model_configuration_file_v1 import ( + 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)) + + +@pytest.fixture +def dataset(tmp_path): + return Dataset(tmp_path) + + +@pytest.fixture(scope="class") +def calibrated(tmp_path_factory): + """One calibration, shared by the assertions that inspect its artifacts.""" + data = Dataset(tmp_path_factory.mktemp("calibration")) + return data, data.calibrate() + + +class TestCalibration: + @pytest.mark.unit + def test_the_configuration_that_produced_the_observations_is_the_optimum(self, calibrated): + data, result = calibrated + + # The observed series is what the configuration itself wrote, so the + # global optimum of the search is the configuration, with NSE 1. + assert result.best_objective == pytest.approx(0.0, abs=1e-9) + assert result.best_nse == pytest.approx(1.0, abs=1e-12) + assert result.best_parameters == pytest.approx(data.parameters, abs=1e-12) + assert set(result.best_parameters) == set(CALIBRATION_PARAMETERS) + + @pytest.mark.unit + def test_every_evaluation_is_a_row_of_the_table(self, calibrated): + data, result = calibrated + rows = read_evaluations(result.evaluations_csv) + + assert list(rows[0]) == list(EVALUATION_COLUMNS) + assert len(rows) == result.evaluations + assert result.best_objective == pytest.approx(min(float(row["objective"]) for row in rows)) + assert {row["id"] for row in rows} == {row["id"] for row in rows if row["id"]} + assert len(rows) >= len(INIT) + # The linear constraint keeps the candidates whose weights do not add up + # away from the workers, so no evaluation is ever rejected as + # inadmissible and none of them failed. + assert [row["error"] for row in rows if row["error"]] == [] + + @pytest.mark.unit + def test_the_evaluations_ran_in_other_processes_and_none_lingers(self, calibrated): + _, result = calibrated + rows = read_evaluations(result.evaluations_csv) + + pids = [int(row["pid"]) for row in rows] + + assert pids, "at least one evaluation was recorded" + assert os.getpid() not in pids + # One task per child: PCRaster keeps its clone and its rasters for the + # life of the process, so no two evaluations may share one. + assert len(set(pids)) == len(pids) + assert multiprocessing.active_children() == [] + + @pytest.mark.unit + def test_the_summary_carries_the_settings_and_the_population_size(self, calibrated): + data, result = calibrated + summary = json.loads(result.result_json.read_text(encoding="utf-8")) + + assert summary["nfev"] == result.evaluations + assert summary["nit"] == result.generations + assert summary["population_size"] == len(INIT) + assert summary["best_objective"] == pytest.approx(result.best_objective) + assert summary["best_nse"] == pytest.approx(result.best_nse) + assert summary["best_parameters"] == pytest.approx(result.best_parameters) + assert summary["settings"]["seed"] == 1 + assert summary["settings"]["workers"] == 2 + assert summary["settings"]["variable"] == "arn" + assert summary["settings"]["maxiter"] == 1 + assert summary["settings"]["polish"] is False + assert summary["settings"]["temp_dir"] == str(data.temp_dir) + assert summary["settings"]["init"] == INIT.tolist() + + @pytest.mark.unit + def test_the_calibrated_configuration_loads_with_the_best_parameters(self, calibrated): + data, result = calibrated + + assert result.calibrated_config.name == "config-calibrated.json" + configuration = ModelConfiguration(result.calibrated_config, validate_input=False) + assert configuration.calibration_parameters.model_dump() == pytest.approx( + result.best_parameters, abs=1e-12 + ) + document = json.loads(result.calibrated_config.read_text(encoding="utf-8")) + assert "CALIBRATION" in document, "a legacy configuration stays legacy" + + @pytest.mark.unit + def test_the_best_row_of_the_table_is_what_was_written_back(self, calibrated): + _, result = calibrated + rows = read_evaluations(result.evaluations_csv) + best = min(rows, key=lambda row: float(row["objective"])) + parameters = {name: float(best[name]) for name in CALIBRATION_PARAMETERS} + + # The optimum SciPy reports and the best record of the run are one and + # the same evaluation, down to the last bit, and it is that evaluation + # the calibrated configuration carries. + assert float(best["objective"]) == result.best_objective + assert float(best["nse"]) == result.best_nse + assert parameters == result.best_parameters + configuration = ModelConfiguration(result.calibrated_config, validate_input=False) + assert configuration.calibration_parameters.model_dump() == pytest.approx( + parameters, abs=1e-12 + ) + + @pytest.mark.unit + def test_the_temporary_directory_is_left_empty(self, calibrated): + data, _ = calibrated + + assert list(data.temp_dir.iterdir()) == [] + + @pytest.mark.unit + def test_the_run_directory_holds_the_three_artifacts(self, calibrated): + _, result = calibrated + + assert result.run_dir.is_dir() + assert result.evaluations_csv.is_file() + assert result.result_json.is_file() + assert result.calibrated_config.is_file() + assert result.message + assert result.generations == 1 + + +class _RecordingExecutor(ProcessPoolExecutor): + """A process pool that remembers how it was built and when it was closed.""" + + built: list["_RecordingExecutor"] = [] + + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.built_with = kwargs + self.shutdowns = 0 + _RecordingExecutor.built.append(self) + + def shutdown(self, *args, **kwargs): + self.shutdowns += 1 + super().shutdown(*args, **kwargs) + + +@pytest.fixture +def broken_search(monkeypatch): + """Record what the search is given, and let it fail without running a model.""" + call = {} + + def explode(function, search_bounds, **kwargs): + call["function"] = function + call["bounds"] = search_bounds + call.update(kwargs) + raise RuntimeError("the search broke") + + _RecordingExecutor.built.clear() + monkeypatch.setattr(scipy_optimize, "differential_evolution", explode) + monkeypatch.setattr(runner, "ProcessPoolExecutor", _RecordingExecutor) + return call + + +class TestProcessTopology: + @pytest.mark.unit + def test_the_pool_is_a_spawned_one_and_is_closed_when_the_search_fails( + self, dataset, broken_search + ): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate() + + (executor,) = _RecordingExecutor.built + assert executor.built_with["max_workers"] == 2 + assert executor.built_with["max_tasks_per_child"] == 1 + assert executor.built_with["mp_context"].get_start_method() == "spawn" + # The pool is closed even though the search raised, and the map of that + # very pool is what the search was told to evaluate with. + assert executor.shutdowns == 1 + assert broken_search["workers"].__self__ is executor + assert multiprocessing.active_children() == [] + + @pytest.mark.unit + def test_the_search_is_given_the_documented_arguments(self, dataset, broken_search): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate() + + assert broken_search["bounds"] == bounds() + assert broken_search["strategy"] == "best1exp" + assert broken_search["maxiter"] == 1 + assert broken_search["popsize"] == 5 + assert broken_search["mutation"] == (0.5, 1.0) + assert broken_search["recombination"] == 0.7 + assert broken_search["updating"] == "deferred" + assert broken_search["polish"] is False + assert broken_search["rng"] == 1 + assert broken_search["init"] is INIT + # The configuration under calibration is the first member of the + # population: SciPy puts x0 in place of the first row of init. + assert broken_search["x0"] == pytest.approx(parameters_to_vector(dataset.parameters)) + + @pytest.mark.unit + def test_the_search_is_constrained_to_weights_that_add_up_to_at_most_one( + self, dataset, broken_search + ): + with pytest.raises(RuntimeError, match="the search broke"): + dataset.calibrate() + + constraint = broken_search["constraints"] + coefficients = np.atleast_2d(constraint.A) + expected = [ + [1.0 if name in ("w_1", "w_2") else 0.0 for name in FREE_PARAMETERS], + ] + + assert coefficients.tolist() == expected + assert np.atleast_1d(constraint.ub).tolist() == [1.0] + assert np.all(np.isneginf(np.atleast_1d(constraint.lb))) + + +class TestProgress: + @pytest.mark.unit + def test_the_callback_takes_the_argument_scipy_passes_by_keyword(self): + # SciPy reads the signature of the callback: only a callback whose one + # parameter is called ``intermediate_result`` is given the result of the + # generation, any other is called with the old ``(x, convergence)`` pair. + assert set(inspect.signature(_Progress(Path())).parameters) == {"intermediate_result"} + + @pytest.mark.unit + def test_the_callback_logs_the_generation_and_lets_the_search_go_on( + self, tmp_path, caplog, capfd + ): + evaluations = tmp_path / "evaluations" + evaluations.mkdir() + (evaluations / "one.json").write_text("{}", encoding="utf8") + progress = _Progress(evaluations) + + with caplog.at_level(logging.INFO, logger="rubem.calibration.runner"): + stop = progress(scipy_optimize.OptimizeResult(fun=12.5, nit=3)) + + assert stop is False, "the callback never halts the search" + assert progress.generations == 3 + assert "Generation 3" in caplog.text + assert "12.5" in caplog.text + assert "1 evaluation" in caplog.text + assert capfd.readouterr().out == "", "the library logs, it does not print" + + +class TestReproducibility: + @pytest.mark.unit + def test_two_runs_of_one_seed_write_the_same_table(self, dataset, capfd): + first = calibrate( + dataset.config_file, + dataset.observed, + dataset.run_dir / "first", + dataset.settings(), + ) + second = calibrate( + dataset.config_file, + dataset.observed, + dataset.run_dir / "second", + dataset.settings(), + ) + + def candidates(result): + return [ + tuple(row[name] for name in CALIBRATION_PARAMETERS) + for row in read_evaluations(result.evaluations_csv) + ] + + # Same seed, same initial population, same candidates, in the same + # order: the 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_the_stations_are_the_values_of_the_sample_locations_raster(self, dataset): + samples = ModelConfiguration( + dataset.config, validate_input=False + ).raster_files.sample_locations + + # The no-data cells of the nominal map are not stations, and the ids are + # the strings the model writes in the header of its table. + assert _station_ids(samples) == {"1", "2"} + + @pytest.mark.unit + def test_a_background_of_zeroes_is_not_a_station(self, tmp_path): + from osgeo import gdal + + gdal.UseExceptions() + gdal.AllRegister() + raster = tmp_path / "samples-with-zeroes.tif" + handle = gdal.GetDriverByName("GTiff").Create(str(raster), 3, 3, 1, gdal.GDT_Int32) + try: + values = np.zeros((3, 3), dtype=np.int32) + values[0, 0], values[2, 2] = 1, 2 + handle.GetRasterBand(1).WriteArray(values) + handle.FlushCache() + finally: + handle = None + + # A raster may mark its background with zeroes and declare no no-data + # value at all; the validation of the inputs counts neither the zeroes + # nor the missing cells as a sample, and neither does this. + assert _station_ids(raster) == {"1", "2"} + + @pytest.mark.unit + def test_the_stations_of_a_zones_aggregation_are_not_checked(self): + # The model remaps the zone ids to 1..N at run time, so the parent has + # no station id to compare and the check only covers the steps. + observed = Series(steps=np.array([1, 2, 3]), stations={"A": np.zeros(3)}) + + assert ( + _check_observed_series( + observed, + first_step=1, + last_step=3, + spinup_steps=0, + aggregation=Aggregation.ZONES, + station_ids={"1", "2"}, + ) + is None + ) + + +class TestSettings: + @pytest.mark.unit + def test_the_default_decision_vector_is_the_configuration_under_calibration(self, dataset): + vector = parameters_to_vector(dataset.parameters) + + assert vector == pytest.approx([4.5, 0.5, 0.333, 0.333, 5.0, 0.5, 0.5, 0.5]) + + @pytest.mark.unit + def test_the_defaults_are_the_documented_ones(self): + settings = CalibrationSettings() + + assert (settings.variable, settings.spinup_steps) == ("arn", 0) + assert (settings.maxiter, settings.popsize) == (100, 15) + assert settings.init == "sobol" + assert settings.strategy == "best1exp" + assert settings.mutation == (0.5, 1.0) + assert settings.recombination == 0.7 + assert settings.polish is False + assert settings.workers is None + + +@pytest.fixture(scope="class") +def calibrated_v1(tmp_path_factory): + """One calibration of the same dataset written as a format 1.0 configuration.""" + tmp_path = tmp_path_factory.mktemp("calibration-v1") + data = Dataset(tmp_path) + legacy = ModelConfigurationFile.model_validate(data.config) + data.config_file = tmp_path / "basin.json" + data.config_file.write_text( + json.dumps(ModelConfigurationFileV1.from_legacy(legacy).to_dict()), encoding="utf8" + ) + return data, data.calibrate() + + +class TestFormatV1Configuration: + @pytest.mark.unit + def test_a_format_1_0_configuration_is_calibrated_and_written_back_in_its_format( + self, calibrated_v1 + ): + data, result = calibrated_v1 + + assert result.best_objective == pytest.approx(0.0, abs=1e-9) + assert result.best_parameters == pytest.approx(data.parameters, abs=1e-12) + assert result.calibrated_config.name == "basin-calibrated.json" + document = json.loads(result.calibrated_config.read_text(encoding="utf-8")) + assert document["version"] == "1.0" + assert document["model_calibration_parameters"]["b"] == pytest.approx( + result.best_parameters["beta"], abs=1e-12 + ) + configuration = ModelConfiguration(result.calibrated_config, validate_input=False) + assert configuration.calibration_parameters.model_dump() == pytest.approx( + result.best_parameters, abs=1e-12 + ) + + +class TestEvaluationBudget: + @pytest.mark.unit + def test_a_sobol_population_is_rounded_up_to_a_power_of_two(self): + # The documented arithmetic: max(5, popsize * 8) members, rounded up to + # the next power of two, because a Sobol' sequence is balanced only over + # a power-of-two sample. + assert _population_size(CalibrationSettings()) == 128 + assert _population_size(CalibrationSettings(popsize=1)) == 8 + assert _population_size(CalibrationSettings(popsize=16)) == 128 + + @pytest.mark.unit + def test_another_initialization_uses_the_members_it_asks_for(self): + assert _population_size(CalibrationSettings(init="latinhypercube", popsize=3)) == 24 + assert _population_size(CalibrationSettings(init="random", popsize=15)) == 120 + assert _population_size(CalibrationSettings(init=INIT)) == len(INIT) + + @pytest.mark.unit + @pytest.mark.parametrize("popsize", [1, 5, 15, 16]) + @pytest.mark.parametrize("init", ["sobol", "latinhypercube", "random"]) + def test_the_announced_population_is_the_one_scipy_builds(self, popsize, init): + # The announced number is what the command line reports as the budget, + # so it has to be SciPy's own arithmetic and not an approximation of it. + optimum = scipy_optimize.differential_evolution( + lambda vector: float(np.sum(vector)), + bounds(), + init=init, + popsize=popsize, + maxiter=0, + polish=False, + rng=0, + ) + + assert len(optimum.population) == _population_size( + CalibrationSettings(popsize=popsize, init=init) + ) + + @pytest.mark.unit + def test_the_summary_reports_the_population_the_search_really_built(self, dataset, monkeypatch): + # A named initialization resizes the population, so the summary must + # report the population of the result and not the number this module + # announced for the budget. + def stub(function, search_bounds, **kwargs): + population = np.asarray(kwargs["init"], dtype=np.float64) + energies = np.asarray([function(member) for member in population[:2]]) + enlarged = np.repeat(population, 2, axis=0) + return scipy_optimize.OptimizeResult( + x=population[0], + fun=float(energies[0]), + nfev=len(energies), + nit=1, + success=True, + message="stub", + population=enlarged, + population_energies=np.full(len(enlarged), np.inf), + ) + + monkeypatch.setattr(scipy_optimize, "differential_evolution", stub) + + result = dataset.calibrate() + summary = json.loads(result.result_json.read_text(encoding="utf-8")) + + assert summary["population_size"] == 2 * len(INIT) + assert _population_size(dataset.settings()) == len(INIT) + + +class TestOptionalDependency: + @pytest.mark.unit + def test_importing_the_package_does_not_import_scipy(self, tmp_path): + # An installation without the rubem[calibration] extra must still be + # able to import the package, so SciPy is imported inside the functions + # that need it. A fresh interpreter is the only honest probe: the test + # session has SciPy loaded already. + probe = tmp_path / "probe.py" + probe.write_text( + "import sys\n" + "import rubem.calibration.runner\n" + "import rubem.calibration._worker\n" + "import rubem.calibration.parameters\n" + "print('scipy' in sys.modules)\n", + encoding="utf8", + ) + + environment = {**os.environ, "PYTHONPATH": REPO_ROOT} + completed = subprocess.run( + [sys.executable, str(probe)], + capture_output=True, + text=True, + check=True, + cwd=REPO_ROOT, + env=environment, + ) + + assert completed.stdout.strip() == "False" diff --git a/tests/unit/calibration/test_worker.py b/tests/unit/calibration/test_worker.py new file mode 100644 index 00000000..301d1011 --- /dev/null +++ b/tests/unit/calibration/test_worker.py @@ -0,0 +1,220 @@ +"""One evaluation of a candidate: the model run, the record and the failure modes.""" + +import json +import os +import shutil +from pathlib import Path + +import pytest + +from rubem.api import Model +from rubem.calibration._worker import EvaluationContext, _derived_document, evaluate +from rubem.calibration.objective import INADMISSIBLE_OBJECTIVE, read_series +from rubem.calibration.parameters import ( + FREE_PARAMETERS, + parameters_to_vector, + vector_to_parameters, +) +from rubem.configuration.model_configuration import ModelConfiguration +from rubem.configuration.model_configuration_file_v1 import ( + VARIABLE_IDS, + ModelConfigurationFileV1, +) +from tests.helpers.synthetic import write_synthetic_dataset + +TIMESTEPS = 3 + + +class Dataset: + """The synthetic dataset, its configuration and the series a plain run wrote.""" + + def __init__(self, tmp_path): + self.config = write_synthetic_dataset(str(tmp_path), timesteps=TIMESTEPS) + self.observed_path = Model.from_config(self.config).run().time_series["arn"][0] + self.configuration = ModelConfiguration(self.config, validate_input=False) + self.document = ModelConfigurationFileV1.from_legacy(self.configuration.file).to_dict() + self.vector = parameters_to_vector(self.configuration.calibration_parameters.model_dump()) + self.temp_dir = tmp_path / "temp" + self.evaluations_dir = tmp_path / "evaluations" + self.temp_dir.mkdir() + self.evaluations_dir.mkdir() + + def context(self, observed_path=None, variable="arn", spinup_steps=0): + """An evaluation context reading ``observed_path``, the plain run by default.""" + return EvaluationContext( + document=self.document, + base_dir=self.configuration.base_dir, + variable=variable, + observed=read_series(observed_path or self.observed_path), + spinup_steps=spinup_steps, + temp_dir=str(self.temp_dir), + evaluations_dir=str(self.evaluations_dir), + validate_input=False, + ) + + def records(self): + """Every record written so far, without the files still being staged.""" + return [ + json.loads(path.read_text(encoding="utf-8")) + for path in sorted(self.evaluations_dir.glob("*.json")) + if not path.name.startswith(".") + ] + + def staged(self): + """The staging files left behind by an interrupted write, if any.""" + return [path.name for path in self.evaluations_dir.iterdir() if ".tmp" in path.name] + + +@pytest.fixture +def dataset(tmp_path): + return Dataset(tmp_path) + + +class TestSuccessfulEvaluation: + @pytest.mark.unit + def test_the_configuration_of_the_observations_scores_a_perfect_objective(self, dataset): + # The observed series is what the configuration itself produced, so its + # own parameters reproduce it exactly: NSE 1, objective 0. + value = evaluate(dataset.vector, dataset.context()) + + assert value == pytest.approx(0.0, abs=1e-9) + (record,) = dataset.records() + assert record["error"] is None + assert record["nse"] == pytest.approx(1.0, abs=1e-12) + assert record["objective"] == pytest.approx(0.0, abs=1e-9) + assert record["station_nse"] == {"1": pytest.approx(1.0), "2": pytest.approx(1.0)} + assert record["parameters"]["w_3"] == pytest.approx(0.334, abs=1e-12) + assert record["pid"] == os.getpid() + assert record["elapsed_seconds"] >= 0.0 + + @pytest.mark.unit + def test_the_temporary_output_directory_is_removed_after_the_run(self, dataset): + evaluate(dataset.vector, dataset.context()) + + assert list(dataset.temp_dir.iterdir()) == [] + + @pytest.mark.unit + def test_the_record_is_installed_by_a_rename_and_leaves_no_staging_file(self, dataset): + evaluate(dataset.vector, dataset.context()) + + assert dataset.staged() == [] + assert len(dataset.records()) == 1 + + @pytest.mark.unit + def test_a_different_candidate_is_evaluated_and_scores_worse(self, dataset): + vector = dataset.vector.copy() + vector[FREE_PARAMETERS.index("alpha_gw")] = 0.01 + + value = evaluate(vector, dataset.context()) + + assert 0.0 < value < INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert record["error"] is None + assert record["nse"] < 1.0 + + +class TestDerivedConfiguration: + @pytest.mark.unit + def test_it_disables_every_raster_series_and_keeps_one_time_series(self, dataset): + document = _derived_document( + dataset.context(), vector_to_parameters(dataset.vector), "/output" + ) + output = document["model_simulation_output"] + rasters = output["raster_series"] + samples = output["time_series_samples"] + + assert output["dir_path"] == "/output" + assert rasters["formats"] == [] + assert [name for name in VARIABLE_IDS if rasters.get(name)] == [] + assert [name for name in VARIABLE_IDS if samples.get(name)] == ["arn"] + assert samples["formats"] == ["CSV"] + assert samples["aggregation"] == "point" + # What the worker hands to the model has to be a configuration again. + ModelConfigurationFileV1.model_validate(document) + + @pytest.mark.unit + def test_the_run_writes_nothing_but_the_table_of_the_calibrated_variable( + self, dataset, monkeypatch + ): + # The output directory is listed on its way out, so that the listing is + # of the directory the evaluation really wrote and the removal still + # happens exactly as it does without the test. + listed = [] + removal = shutil.rmtree + + def spy(path, **kwargs): + listed.extend(sorted(entry.name for entry in Path(path).rglob("*"))) + removal(path, **kwargs) + + monkeypatch.setattr(shutil, "rmtree", spy) + + evaluate(dataset.vector, dataset.context()) + + assert listed == ["metadata.json", "tss_arn.csv"] + + +class TestRejectedCandidate: + @pytest.mark.unit + def test_an_inadmissible_candidate_is_recorded_without_a_model_run(self, dataset): + vector = dataset.vector.copy() + vector[FREE_PARAMETERS.index("w_1")] = 0.6 + vector[FREE_PARAMETERS.index("w_2")] = 0.6 + + value = evaluate(vector, dataset.context()) + + assert value == INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert record["error"] == "inadmissible" + assert record["nse"] is None + assert record["objective"] == INADMISSIBLE_OBJECTIVE + assert list(dataset.temp_dir.iterdir()) == [], "no output directory was created" + + +class TestRecord: + @pytest.mark.unit + def test_the_record_is_staged_in_the_directory_it_is_renamed_into(self, dataset, monkeypatch): + # A rename is atomic only inside one filesystem, so the staging file has + # to be a sibling of the record and not a file in the system temporary + # directory, which may be a different mount. + renames = [] + rename = Path.replace + + def spy(self, target): + renames.append((Path(self), Path(target))) + return rename(self, target) + + monkeypatch.setattr(Path, "replace", spy) + + evaluate(dataset.vector, dataset.context()) + + installed = [pair for pair in renames if pair[1].parent == dataset.evaluations_dir] + assert len(installed) == 1, "the record was installed by exactly one rename" + temporary, target = installed[0] + assert temporary.parent == dataset.evaluations_dir + assert temporary.name.startswith(".") + assert target.name.endswith(".json") + + +class TestFailedEvaluation: + @pytest.mark.unit + def test_observations_of_another_station_are_recorded_as_the_error(self, dataset, tmp_path): + foreign = tmp_path / "foreign.csv" + foreign.write_text("0;A;B\n1;1.0;2.0\n2;2.0;3.0\n3;3.0;4.0\n", encoding="utf8") + + value = evaluate(dataset.vector, dataset.context(observed_path=foreign)) + + assert value == INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert "no station in common" in record["error"] + assert record["error"].startswith("ValueError: ") + assert record["nse"] is None + assert list(dataset.temp_dir.iterdir()) == [], "the output directory was removed anyway" + + @pytest.mark.unit + def test_a_variable_the_run_does_not_write_is_recorded_as_the_error(self, dataset): + value = evaluate(dataset.vector, dataset.context(variable="nope")) + + assert value == INADMISSIBLE_OBJECTIVE + (record,) = dataset.records() + assert record["error"] + assert list(dataset.temp_dir.iterdir()) == [] diff --git a/tests/unit/test_cli_calibrate.py b/tests/unit/test_cli_calibrate.py new file mode 100644 index 00000000..3fb9f4f0 --- /dev/null +++ b/tests/unit/test_cli_calibrate.py @@ -0,0 +1,479 @@ +import json +from pathlib import Path + +import pytest + +from rubem.cli import main +from tests.helpers.synthetic import write_synthetic_dataset + + +@pytest.fixture(name="config_path") +def config_path_fixture(tmp_path): + path = tmp_path / "config.json" + path.write_text(json.dumps(write_synthetic_dataset(str(tmp_path))), encoding="utf8") + return path + + +@pytest.fixture(name="observed_path") +def observed_path_fixture(tmp_path): + path = tmp_path / "observed.csv" + path.write_text("0;1\n1;2.0\n2;3.0\n", encoding="utf8") + return path + + +class FakeCalibration: + """Stand in for ``rubem.calibration.runner.calibrate``. + + It records the arguments the command line built and answers with a result + whose files do not have to exist: the command only prints their paths. + """ + + def __init__(self, run_dir, best_nse=0.87): + from rubem.calibration.runner import CalibrationResult + + self.calls = [] + self.result = CalibrationResult( + best_parameters={ + "alpha": 4.5, + "beta": 0.5, + "w_1": 0.333, + "w_2": 0.333, + "w_3": 0.334, + "rcd": 5.0, + "f": 0.5, + "alpha_gw": 0.5, + "x": 0.5, + }, + best_nse=best_nse, + best_objective=169000.0, + evaluations=128, + generations=2, + success=True, + message="Optimization terminated successfully.", + run_dir=Path(run_dir), + evaluations_csv=Path(run_dir) / "evaluations.csv", + result_json=Path(run_dir) / "result.json", + calibrated_config=Path(run_dir) / "config-calibrated.json", + ) + + def __call__(self, config_path, observed_path, run_dir, settings): + self.calls.append((config_path, observed_path, run_dir, settings)) + return self.result + + +@pytest.fixture(name="fake_calibration") +def fake_calibration_fixture(tmp_path, monkeypatch): + fake = FakeCalibration(tmp_path / "calibration") + monkeypatch.setattr("rubem.calibration.runner.calibrate", fake) + return fake + + +class TestCliCalibrateArgumentsReachTheRunner: + @pytest.mark.unit + def test_the_paths_and_the_documented_defaults( + self, config_path, observed_path, tmp_path, fake_calibration, restore_logging + ): + run_dir = tmp_path / "calibration" + + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(run_dir), + ] + ) + + (given_config, given_observed, given_run_dir, settings) = fake_calibration.calls[0] + assert Path(given_config) == config_path + assert Path(given_observed) == observed_path + assert Path(given_run_dir) == run_dir + assert settings.variable == "arn" + assert settings.spinup_steps == 0 + assert settings.maxiter == 100 + assert settings.popsize == 15 + assert settings.seed is None + assert settings.temp_dir is None + + @pytest.mark.unit + def test_every_option_reaches_the_settings( + self, config_path, observed_path, tmp_path, fake_calibration, restore_logging + ): + temp_dir = tmp_path / "scratch" + + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + "--variable", + "rnf", + "--spinup-steps", + "12", + "--maxiter", + "3", + "--popsize", + "4", + "--seed", + "7", + "--workers", + "2", + "--temp-dir", + str(temp_dir), + ] + ) + + settings = fake_calibration.calls[0][3] + assert settings.variable == "rnf" + assert settings.spinup_steps == 12 + assert settings.maxiter == 3 + assert settings.popsize == 4 + assert settings.seed == 7 + assert settings.workers == 2 + assert Path(settings.temp_dir) == temp_dir + + @pytest.mark.unit + def test_the_runner_keeps_its_own_default_worker_count( + self, config_path, observed_path, tmp_path, fake_calibration, monkeypatch, restore_logging + ): + """Without ``--workers`` the command must not override the runner's default. + + The settings class is replaced by one whose default differs from the + ``None`` the option carries, so forwarding the option unconditionally + would be visible here. + """ + from dataclasses import dataclass + + from rubem.calibration.runner import CalibrationSettings + + @dataclass(frozen=True) + class SettingsWithAWorkerDefault(CalibrationSettings): + workers: int | None = 7 + + monkeypatch.setattr( + "rubem.calibration.runner.CalibrationSettings", SettingsWithAWorkerDefault + ) + + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + assert fake_calibration.calls[0][3].workers == 7 + + +class TestCliCalibrateOutput: + @pytest.mark.unit + def test_the_documented_progress_and_summary_are_printed( + self, config_path, observed_path, tmp_path, fake_calibration, capsys, restore_logging + ): + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + output = capsys.readouterr().out + result = fake_calibration.result + for expected in ( + "Loading configuration and validating inputs...", + "Calibration started...", + "Calibration finished successfully!", + "Best NSE: 0.87", + "Best objective: 169000", + "Evaluations: 128", + "Generations: 2", + str(result.evaluations_csv), + str(result.result_json), + str(result.calibrated_config), + ): + assert expected in output, f"missing line: {expected!r}\n{output}" + + @pytest.mark.unit + def test_a_run_without_a_usable_nse_says_so( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + fake = FakeCalibration(tmp_path / "calibration", best_nse=None) + monkeypatch.setattr("rubem.calibration.runner.calibrate", fake) + + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + assert "Best NSE: n/a" in capsys.readouterr().out + + +class TestCliCalibrateFailures: + @pytest.mark.unit + def test_a_calibration_error_exits_with_one_without_a_traceback( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + from rubem.calibration.runner import CalibrationError + + def fail(*args, **kwargs): + raise CalibrationError("every evaluation failed: no such variable") + + monkeypatch.setattr("rubem.calibration.runner.calibrate", fail) + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + captured = capsys.readouterr().err + assert error.value.code == 1 + assert "every evaluation failed: no such variable" in captured + assert "Traceback" not in captured + + @pytest.mark.unit + def test_an_unexpected_failure_logs_one_traceback( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + def fail(*args, **kwargs): + raise RuntimeError("the pool died") + + monkeypatch.setattr("rubem.calibration.runner.calibrate", fail) + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + captured = capsys.readouterr().err + assert error.value.code == 1 + assert captured.count("Traceback (most recent call last)") == 1 + assert "RUBEM unexpectedly quit." in captured + + @pytest.mark.unit + def test_a_missing_scipy_exits_with_the_installation_guidance( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + monkeypatch.setattr("rubem._deps.missing_calibration_deps", lambda: ["scipy"]) + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + message = str(error.value) + assert "scipy" in message + assert 'pip install "rubem[calibration]"' in message + + @pytest.mark.unit + def test_an_interrupted_calibration_exits_with_two( + self, config_path, observed_path, tmp_path, monkeypatch, capsys, restore_logging + ): + def interrupt(*args, **kwargs): + raise KeyboardInterrupt + + monkeypatch.setattr("rubem.calibration.runner.calibrate", interrupt) + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + ] + ) + + captured = capsys.readouterr().err + assert error.value.code == 2 + assert "RUBEM was interrupted by the user." in captured + assert "Traceback" not in captured + + @pytest.mark.unit + @pytest.mark.parametrize( + ("option", "value"), + [ + ("--maxiter", "0"), + ("--maxiter", "-1"), + ("--popsize", "0"), + ("--popsize", "-1"), + ("--workers", "0"), + ("--workers", "-1"), + ("--spinup-steps", "-1"), + ], + ) + def test_a_search_size_outside_its_range_is_rejected( + self, config_path, observed_path, tmp_path, fake_calibration, option, value, restore_logging + ): + """A generation count, a population and a worker count all start at one.""" + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + option, + value, + ] + ) + + assert error.value.code == 2 + assert not fake_calibration.calls + + @pytest.mark.unit + def test_a_negative_spin_up_is_rejected( + self, config_path, observed_path, tmp_path, fake_calibration, capsys, restore_logging + ): + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(observed_path), + "-o", + str(tmp_path / "calibration"), + "--spinup-steps", + "-1", + ] + ) + + assert error.value.code == 2 + assert not fake_calibration.calls + + @pytest.mark.unit + def test_a_missing_observed_file_is_rejected_before_anything_runs( + self, config_path, tmp_path, fake_calibration, capsys, restore_logging + ): + missing = tmp_path / "missing.csv" + + with pytest.raises(SystemExit) as error: + main( + [ + "calibrate", + "-c", + str(config_path), + "--observed", + str(missing), + "-o", + str(tmp_path / "calibration"), + ] + ) + + assert error.value.code == 2 + assert not fake_calibration.calls + assert "Loading configuration" not in capsys.readouterr().out + + @pytest.mark.unit + def test_an_incomplete_command_is_not_mapped_onto_the_legacy_spelling( + self, config_path, capsys, restore_logging + ): + """``calibrate`` is a subcommand, not the deprecated bare ``-c`` form.""" + with pytest.raises(SystemExit) as error: + main(["calibrate", "-c", str(config_path)]) + + captured = capsys.readouterr() + assert error.value.code == 2 + assert "Missing option '--observed'" in captured.err + assert "is deprecated" not in captured.err + + +class TestCliCalibrateHelp: + @pytest.mark.unit + def test_the_command_is_listed_in_the_help(self, capsys, restore_logging): + with pytest.raises(SystemExit) as error: + main(["--help"]) + + assert error.value.code == 0 + assert "calibrate" in capsys.readouterr().out + + @pytest.mark.unit + def test_the_help_documents_every_option(self, capsys, restore_logging): + with pytest.raises(SystemExit) as error: + main(["calibrate", "--help"]) + + assert error.value.code == 0 + output = capsys.readouterr().out + assert "Calibrate the model parameters against an observed series." in output + for option in ( + "--configfile", + "--observed", + "--run-dir", + "--variable", + "--spinup-steps", + "--maxiter", + "--popsize", + "--seed", + "--workers", + "--temp-dir", + ): + assert option in output, f"missing option: {option}\n{output}" + + @pytest.mark.unit + def test_importing_the_command_line_does_not_need_scipy(self): + """The CLI module stays importable without the optional extra.""" + import subprocess + import sys + + completed = subprocess.run( + [sys.executable, "-c", "import sys, rubem.cli; print('scipy' in sys.modules)"], + cwd=Path(__file__).resolve().parents[2], + capture_output=True, + text=True, + check=True, + ) + assert completed.stdout.strip() == "False" diff --git a/tests/unit/test_deps.py b/tests/unit/test_deps.py index 410d5bca..41bcb5ae 100644 --- a/tests/unit/test_deps.py +++ b/tests/unit/test_deps.py @@ -18,3 +18,54 @@ def test_missing_dependencies_are_named_with_an_actionable_source(self, monkeypa assert name in message assert "conda-forge" in message assert _deps._ENVIRONMENT_YML_URL in message + + +class TestRequireCalibrationDeps: + @pytest.mark.unit + def test_present_dependencies_do_not_raise(self): + _deps.require_calibration_deps() + + @pytest.mark.unit + def test_nothing_is_missing_when_every_dependency_is_importable(self): + assert _deps.missing_calibration_deps() == [] + + @pytest.mark.unit + def test_an_unimportable_dependency_is_reported_as_missing(self, monkeypatch): + monkeypatch.setattr(_deps.importlib.util, "find_spec", lambda name: None) + + assert _deps.missing_calibration_deps() == list(_deps._CALIBRATION_DEPENDENCIES) + + @pytest.mark.unit + def test_a_broken_dependency_counts_as_missing(self, monkeypatch): + def raise_import_error(name): + raise ImportError(name) + + monkeypatch.setattr(_deps.importlib.util, "find_spec", raise_import_error) + + assert _deps.missing_calibration_deps() == list(_deps._CALIBRATION_DEPENDENCIES) + + @pytest.mark.unit + def test_the_message_names_the_dependency_and_the_extra(self): + message = _deps.calibration_deps_message(["scipy"]) + + assert "scipy" in message + assert 'pip install "rubem[calibration]"' in message + + @pytest.mark.unit + def test_the_message_defaults_to_every_calibration_dependency(self): + """The runner asks for the guidance without having computed the list.""" + assert _deps.calibration_deps_message() == _deps.calibration_deps_message( + list(_deps._CALIBRATION_DEPENDENCIES) + ) + + @pytest.mark.unit + def test_missing_dependencies_are_named_with_an_actionable_extra(self, monkeypatch): + monkeypatch.setattr(_deps.importlib.util, "find_spec", lambda name: None) + + with pytest.raises(SystemExit) as error: + _deps.require_calibration_deps() + + message = str(error.value) + for name in _deps._CALIBRATION_DEPENDENCIES: + assert name in message + assert 'pip install "rubem[calibration]"' in message