Support more of LARSIM's branching features

@BGWKlein

We need to add some functionalities from the following LARSIM options to HydPy:

VERZ: simple branching from streams
ERWV: extended branching from streams
VSPE: branching from reservoirs

First thoughts (only on VERZ and ERWV for now):

VERZ seems like a simplification of ERWV. I don't think we need separate HydPy models to reflect both options.

Basically, ERWV works like HydPy-Exch-Branch-HBV96. However, it has the following additional features:

There can be different rules either for different times of the year (e.g. VONBIS 01.06.9999 05:00 01.07.9999 05:00) or for different concrete periods (e.g. VONBIS 01.06.2010 05:00 01.06.2020 05:00).
One can define one time series for each rule. If time series values are not missing, they define the extraction amount (and override the fixed relationship).

Defining rules for concrete periods is an interesting feature, but it's seldom required and isn't currently on our list. So, we can focus on defining rules for different times of the year (seasons). This is something we usually do with the help of the TOYParameter class.

The code for configuring HydPy-Exch-Branch-HBV96 currently looks like this:

from hydpy.models.exch_branch_hbv96 import *
parameterstep()
delta(-1.0)
minimum(-1.0)
xpoints(0.0, 2.0, 4.0, 6.0)
ypoints(output1=[0.0, 1.0, 2.0, 3.0],
        output2=[0.0, 1.0, 0.0, 0.0],
        output3=[0.0, 0.0, 2.0, 6.0])

Delta and Minimum are very special parameters that we needed for compatibility with HBV96.

HydPy-Exch-Branch-HBV96 allows branching to an arbitrary number of places; LARSIM allows branching to only one place.

Extending all these parameters (or at least XPoints and YPoints with "time of year" (toy) functionality seems suboptimal:

xpoints(toy_06=[0.0, 2.0, 4.0, 6.0], toy_07=[0.0, 2.0, 4.0, 6.0])
ypoints(
    toy_06={"output1"=[0.0, 1.0, 2.0, 3.0], "output2"=[0.0, 1.0, 0.0, 0.0], "output3"=[0.0, 0.0, 2.0, 6.0]},
    toy_07={"output1"=[0.0, 1.0, 2.0, 4.0], "output2"=[0.0, 1.0, 0.0, 0.0], "output3"=[0.0, 0.0, 2.0, 6.0]},
)

Keeping the data for each time of year seems preferable. For example;

from hydpy.models.exch_branch_larsim import Rule
rules(
    toy_06=Rule(
        [0.0, 2.0, 4.0, 6.0], 
        output1=[0.0, 1.0, 2.0, 3.0], 
        output2=[0.0, 1.0, 0.0, 0.0],
        output3=[0.0, 0.0, 2.0, 6.0],
    ),
    toy_07=Rule(
        [0.0, 2.0, 4.0, 6.0], 
        output1=[0.0, 1.0, 2.0, 4.0], 
        output2=[0.0, 1.0, 0.0, 0.0],
        output3=[0.0, 0.0, 2.0, 6.0],
    ),
)

My original gut feeling was to extend HydPy-Branch-HBV96 with the new features. However, because Delta and Minimum are a little special and the API we've designed so far differs significantly, I'm now considering creating a new model, HydPy-Branch-LARSIM (working title).

Regarding the externally defined time series, I would prefer to keep things simple. LARSIM lets you assign different time series to individual rules. I want to delegate such work to preprocessing, meaning users would have to compile multiple series into one. Following the mentioned LARSIM convention, nan values in the time series mean that the currently relevant fixed rule applies.

Two problems:

In HydPy, nan usually indicates data that is unfortunately missing. Here, it may mean the same but could also mean "just use the standard rule here". Hence, an activated checkseries flag would often cause false positives. (Didn't we have a similar problem before?)
Users would have to prepare nan time series, even if they only want to rely on fixed rules. To avoid this, we could create two versions of HydPy-Branch-LARSIM or allow adding a submodel that reads the time series only if desired. (I have a small tendency for the latter.)

cc @BGWKlein @parvumflumen

Other HydPy models interpret, for example, toy_06_01=2.0 and toy_07_01=3.0 so that 2.0 applies only on each June 1 and 3.0 applies only on each July 1. The values for all other days in June are linearly interpolated. I think we should stick to this practice. This would give us flexibility, but also make defining a complete LARSIM-like, discontinuous annual cycle more laborious (assuming a simulation step of one day):

rules(
    toy_05_31_12=Rule(
        [0.0, 2.0, 4.0, 6.0], 
        output1=[0.0, 1.0, 2.0, 4.0], 
        output2=[0.0, 1.0, 0.0, 0.0],
        output3=[0.0, 0.0, 2.0, 6.0],
    ),
    toy_06_01_12=Rule(
        [0.0, 2.0, 4.0, 6.0], 
        output1=[0.0, 1.0, 2.0, 3.0], 
        output2=[0.0, 1.0, 0.0, 0.0],
        output3=[0.0, 0.0, 2.0, 6.0],
    ),
    toy_06_30_12=Rule(
        [0.0, 2.0, 4.0, 6.0], 
        output1=[0.0, 1.0, 2.0, 3.0], 
        output2=[0.0, 1.0, 0.0, 0.0],
        output3=[0.0, 0.0, 2.0, 6.0],
    ),
    toy_07_01_12=Rule(
        [0.0, 2.0, 4.0, 6.0], 
        output1=[0.0, 1.0, 2.0, 4.0], 
        output2=[0.0, 1.0, 0.0, 0.0],
        output3=[0.0, 0.0, 2.0, 6.0],
    ),
)

This could eventually be simplified later, for example by automatically reusing equal neighbouring rules in string representations:

rules(
    toy_05_31_12=Rule(
        [0.0, 2.0, 4.0, 6.0], 
        output1=[0.0, 1.0, 2.0, 4.0], 
        output2=[0.0, 1.0, 0.0, 0.0],
        output3=[0.0, 0.0, 2.0, 6.0],
    ),
    toy_06_01_12=Rule(
        [0.0, 2.0, 4.0, 6.0], 
        output1=[0.0, 1.0, 2.0, 3.0], 
        output2=[0.0, 1.0, 0.0, 0.0],
        output3=[0.0, 0.0, 2.0, 6.0],
    ),
    toy_06_30_12="toy_06_01_12"
    toy_07_01_12="toy_05_31_12",
)

Or maybe we could support additional keywords:

rules(
    from_06=Rule(
        [0.0, 2.0, 4.0, 6.0], 
        output1=[0.0, 1.0, 2.0, 4.0], 
        output2=[0.0, 1.0, 0.0, 0.0],
        output3=[0.0, 0.0, 2.0, 6.0],
    ),
    from_07=Rule(
        [0.0, 2.0, 4.0, 6.0], 
        output1=[0.0, 1.0, 2.0, 3.0], 
        output2=[0.0, 1.0, 0.0, 0.0],
        output3=[0.0, 0.0, 2.0, 6.0],
    ),
)

But these steps don't seem urgent.

In HydPy, nan usually indicates data that is unfortunately missing. Here, it may mean the same but could also mean "just use the standard rule here". Hence, an activated checkseries flag would often cause false positives. (Didn't we have a similar problem before?)

I see two simple solutions:

Use a specific negative value to indicate that HydPy-Exch-Branch-LARSIM should use the currently relevant fixed rule (e.g. -999.9).
Let HydPy-Exch-Branch-LARSIM understand any value below zero this way.
Use numpy.inf (and also -numpy.inf).

I think the third option is the best because we can use numpy.isinf to identify inf (and -inf), and normal negative values could more easily lead to misleading analysis results (for example, when calculating averages, whereas numpy.mean([inf, 0.0] is inf).

I am thinking about two possible ways to refine the mentioned "Rule" idea.

`Rules` as a "special" parameter

Then, no Rules instances would be available during simulations. Rules would collect information from the user only. Corresponding derived parameters would convert this information to data usable during simulation runs.

We already discussed officially introducing "special" parameters, which users can access directly like "normal" control parameters, but which are clearly distinguishable from them. Currently, the inheritance tree is a little misleading, giving the impression that all already available special parameters have all the features of normal control parameters, which is not the case.

I think if we implement Rules from the ground up, it would be a perfect opportunity to introduce something like a SpecialParameter class and derive Rules from it. However, subclassing the existing special parameters could introduce breaking changes, but we could postpone that (it would be a cool feature for HydPy 7).

`Rules` subclassed from SeasonalInterpolator

Searching through our existing base classes, I found SeasonalInterpolator. Basically, it does what we want (interpolate for specific times of the year based on some fixed x-y relationship and then interpolate the results temporarily). However, there is a limitation regarding the number of x-y relationships per time of year. We need one input (x) and an arbitrary number of outputs (y). SeasonalInterpolator itself does not seem to be limited in this regard. However, it uses instances of the base class InterpAlgorithm to represent the x-y relationships. Currently, we have only two concrete, applicable subclasses of it. PPoly allows for linear and different kinds of piecewise polynomial interpolation, but unfortunately only for univariate output (one y value). ANN allows for arbitrary relationships and an arbitrary number of outputs, but it relies on neural network equations, which is too much effort to configure in most situations.

So we would have to implement (at least) two enhancements:

Let PPoly support multivariate output or, more likely, implement a container class that handles multiple PPoly instances (or something similar).
If we interpolate multiple relationships independently via splines, water balance violations become likely. So we would either have to suggest preferring linear interpolation (and eventually introduce some kind of restrictions or warnings). Or we would have to let the container class normalise the sum of all y values to the current x values (or something similar).

With respect to the last point, we must keep in mind that HydPy-Exch-Branch-HBV96 actively allows water violations (see the examples above), while LARSIM likely does not (at least, if I understand the documentation correctly).

Rules subclassed from SeasonalInterpolator

I currently favour this option. I have a prototype of PPolys that handles multiple PPoly instances. Later, one should be able to pass one PPolys instance to a parameter derived from SeasonalInterpolator for each time of year. The following examples show my current ideas:

from scipy.interpolate import PchipInterpolator
from hydpy import Poly, PPoly, PPolys
ppolys = PPolys()
ppolys(
    main=PPolys.REST,
    branch1=PPoly.from_data(xs=[0.0, 10.0, 16.0], ys=[0.0, 0.0, 1.0]),
    branch2=PPoly.from_data(
        xs=[0.0, 10.0, 16.0, 22.0, 28.0], ys=[0.0, 0.0, 2.0, 0.0, 0.0]
    ),
    branch3=PPoly.from_data(
        xs=[0.0, 10.0, 16.0, 22.0, 28.0], ys=[0.0, 0.0, 3.0, 12.0, 12.0]
    ),
    supplier=PPoly.from_data(
        xs=[0.0, 1.0, 2.0], ys=[-2.0, 0.0, 0.0], method=PchipInterpolator
    ),
)

The HydPy-Exch-Branch-LARSIM submodel receives inflow from a normal inlet node and connects to five outlet nodes named main, branch1, branch2, branch3, and supplier. This complexity is likely seldom required (or never), but it gives an impression of all possible features I currently see.

For all branches, interpolation happens piecewise linearly. From 10 m³/s on, all additional discharge is routed into the branches. For branch 1, there is linear extrapolation. The routing to branch 2 has a maximum of 2 m³/s for an inflow of 16 m³/s and falls back to zero for larger inflow values. All remaining excess inflow reaches branch 3, until its maximum capacity of 12 m³/s is reached.

The supplier node works differently. First, it sends water when inflow values are small. Second, this supply decreases nonlinearly from 2 m³/s to zero.

With respect to the three branch nodes, one could configure the release to the main river by defining another x-y-based piecewise linear relationship, without violating the water balance. However, this is error-prone and very hard, or even impossible, when working with spline interpolation, as with the supplier node. Hence, we could say that the main branch should always receive the rest, but this might not always be desired (for example, we might want to "lose" water without explicitly sending it to another node). Hence, we could tell the PPolys instance to always calculate the rest by itself and send it to the main node by assigning the constant PPolys.REST. This would also allow more thorough verification (for example, the sum of all y values must not exceed the corresponding x value). Maybe we would want to add other options later by adding more constants.

Not that the value passed to the supplier node only depends on the inflow to the HydPy-Exch-Branch-LARSIM instance. This could be problematic if the supplier node extracts water from another model instance that in fact cannot provide enough water.

Also note that the automatic string representation of the above example looks even more complicated (but this should be fixable):

    PPolys(
        main=PPolys.REST,
        branch1=PPoly(
            Poly(x0=0.0, cs=(0.0, 0.0)),
            Poly(x0=10.0, cs=(0.0, 0.166667)),
        ),
        branch2=PPoly(
            Poly(x0=0.0, cs=(0.0, 0.0)),
            Poly(x0=10.0, cs=(0.0, 0.333333)),
            Poly(x0=16.0, cs=(2.0, -0.333333)),
            Poly(x0=22.0, cs=(0.0, 0.0)),
        ),
        branch3=PPoly(
            Poly(x0=0.0, cs=(0.0, 0.0)),
            Poly(x0=10.0, cs=(0.0, 0.5)),
            Poly(x0=16.0, cs=(3.0, 1.5)),
            Poly(x0=22.0, cs=(12.0, 0.0)),
        ),
        supplier=PPoly(
            Poly(x0=0.0, cs=(-2.0, 3.0, 0.0, -1.0)),
            Poly(x0=1.0, cs=(0.0,)),
        ),
    )

I could get the first preliminary doctest running with relative ease (still using 1 instead of PPolys.REST and with fixed x-y data compared to the example mentioned above) without bumping into an unexpected problem:

    >>> from scipy.interpolate import PchipInterpolator
    >>> from hydpy import Poly, PPoly, PPolys
    >>> ppolys = PPolys()
    >>> ppolys(
    ...     main=PPoly.from_data(
    ...         xs=[0.0, 10.0, 22.0, 28.0], ys=[0.0, 10.0, 10.0, 15.0]
    ...     ),
    ...     branch1=PPoly.from_data(xs=[0.0, 10.0, 16.0], ys=[0.0, 0.0, 1.0]),
    ...     branch2=PPoly.from_data(
    ...         xs=[0.0, 10.0, 16.0, 22.0, 28.0], ys=[0.0, 0.0, 2.0, 0.0, 0.0]
    ...     ),
    ...     branch3=PPoly.from_data(
    ...         xs=[0.0, 10.0, 16.0, 22.0, 28.0], ys=[0.0, 0.0, 3.0, 10.0, 10.0]
    ...     ),
    ... )

    >>> from hydpy import print_vector
    >>> def check(print_derivatives=False):
    ...     for i in [
    ...         0.0, 2.5, 5.0, 7.5, 10.0, 13.0, 16.0, 19.0, 22.0, 25.0, 28.0, 31.0, 34.0
    ...     ]:
    ...         ppolys.inputs[0] = i
    ...         ppolys.calculate_values()
    ...         assert i == sum(ppolys.outputs), f"{i} vs {sum(ppolys.outputs)}"
    ...         print(i, end=": ")
    ...         if print_derivatives:
    ...             ppolys.calculate_derivatives()
    ...             print_vector(ppolys.output_derivatives)
    ...         else:
    ...             print_vector(ppolys.outputs)

    >>> check()
    0.0: 0.0, 0.0, 0.0, 0.0
    2.5: 2.5, 0.0, 0.0, 0.0
    5.0: 5.0, 0.0, 0.0, 0.0
    7.5: 7.5, 0.0, 0.0, 0.0
    10.0: 10.0, 0.0, 0.0, 0.0
    13.0: 10.0, 0.5, 1.0, 1.5
    16.0: 10.0, 1.0, 2.0, 3.0
    19.0: 10.0, 1.5, 1.0, 6.5
    22.0: 10.0, 2.0, 0.0, 10.0
    25.0: 12.5, 2.5, 0.0, 10.0
    28.0: 15.0, 3.0, 0.0, 10.0
    31.0: 17.5, 3.5, 0.0, 10.0
    34.0: 20.0, 4.0, 0.0, 10.0

    >>> check(print_derivatives=True)
    0.0: 1.0, 0.0, 0.0, 0.0
    2.5: 1.0, 0.0, 0.0, 0.0
    5.0: 1.0, 0.0, 0.0, 0.0
    7.5: 1.0, 0.0, 0.0, 0.0
    10.0: 0.0, 0.166667, 0.333333, 0.5
    13.0: 0.0, 0.166667, 0.333333, 0.5
    16.0: 0.0, 0.166667, -0.333333, 1.166667
    19.0: 0.0, 0.166667, -0.333333, 1.166667
    22.0: 0.833333, 0.166667, 0.0, 0.0
    25.0: 0.833333, 0.166667, 0.0, 0.0
    28.0: 0.833333, 0.166667, 0.0, 0.0
    31.0: 0.833333, 0.166667, 0.0, 0.0
    34.0: 0.833333, 0.166667, 0.0, 0.0

    >>> ppolys(
    ...     main=1,
    ...     branch1=PPoly.from_data(xs=[0.0, 10.0, 16.0], ys=[0.0, 0.0, 1.0]),
    ...     branch2=PPoly.from_data(
    ...         xs=[0.0, 10.0, 16.0, 22.0, 28.0], ys=[0.0, 0.0, 2.0, 0.0, 0.0]
    ...     ),
    ...     branch3=PPoly.from_data(
    ...         xs=[0.0, 10.0, 16.0, 22.0, 28.0], ys=[0.0, 0.0, 3.0, 10.0, 10.0]
    ...     ),
    ...     supplier=PPoly.from_data(
    ...         xs=[0.0, 8.0, 10.0], ys=[-2.0, 0.0, 0.0], method=PchipInterpolator
    ...     ),
    ... )

    >>> check()
    0.0: 2.0, 0.0, 0.0, 0.0, -2.0
    2.5: 3.504395, 0.0, 0.0, 0.0, -1.004395
    5.0: 5.316406, 0.0, 0.0, 0.0, -0.316406
    7.5: 7.509277, 0.0, 0.0, 0.0, -0.009277
    10.0: 10.0, 0.0, 0.0, 0.0, 0.0
    13.0: 10.0, 0.5, 1.0, 1.5, 0.0
    16.0: 10.0, 1.0, 2.0, 3.0, 0.0
    19.0: 10.0, 1.5, 1.0, 6.5, 0.0
    22.0: 10.0, 2.0, 0.0, 10.0, 0.0
    25.0: 12.5, 2.5, 0.0, 10.0, 0.0
    28.0: 15.0, 3.0, 0.0, 10.0, 0.0
    31.0: 17.5, 3.5, 0.0, 10.0, 0.0
    34.0: 20.0, 4.0, 0.0, 10.0, 0.0

    >>> check(print_derivatives=True)
    0.0: 0.55, 0.0, 0.0, 0.0, 0.45
    2.5: 0.658398, 0.0, 0.0, 0.0, 0.341602
    5.0: 0.796094, 0.0, 0.0, 0.0, 0.203906
    7.5: 0.963086, 0.0, 0.0, 0.0, 0.036914
    10.0: 0.0, 0.166667, 0.333333, 0.5, 0.0
    13.0: 0.0, 0.166667, 0.333333, 0.5, 0.0
    16.0: 0.0, 0.166667, -0.333333, 1.166667, 0.0
    19.0: 0.0, 0.166667, -0.333333, 1.166667, 0.0
    22.0: 0.833333, 0.166667, 0.0, 0.0, 0.0
    25.0: 0.833333, 0.166667, 0.0, 0.0, 0.0
    28.0: 0.833333, 0.166667, 0.0, 0.0, 0.0
    31.0: 0.833333, 0.166667, 0.0, 0.0, 0.0
    34.0: 0.833333, 0.166667, 0.0, 0.0, 0.0

Many things need some love (for example, the error messages raised by verify), but nothing seems to speak against going this way.

I tried to put everything together quickly and dirty, and in principle, everything seems to work.

The following configuration reproduces the single integration test of HydPy-Exch-Branch-HBV96:

>>> rules(
...     PPolys(
...         output1=PPoly.from_data(xs=[0.0, 2.0], ys=[0.0, 1.0]),
...         output2=PPoly.from_data(xs=[0.0, 2.0, 4.0, 6.0], ys=[0.0, 1.0, 0.0, 0.0]),
...         output3=PPoly.from_data(xs=[0.0, 2.0, 4.0, 6.0], ys=[0.0, 0.0, 2.0, 6.0]),
...     )
... )

    >>> test("exch_branch_larsim_ex1")
    |  date | originalinput |             outputs | input1 | input2 | output1 | output2 | output3 |
    -----------------------------------------------------------------------------------------------
    | 00:00 |          -1.0 | -0.5  -0.5      0.0 |   -2.0 |    1.0 |    -0.5 |    -0.5 |     0.0 |
    | 01:00 |          -1.0 | -0.5  -0.5      0.0 |   -1.0 |    0.0 |    -0.5 |    -0.5 |     0.0 |
    | 02:00 |           0.0 |  0.0   0.0      0.0 |    0.0 |    0.0 |     0.0 |     0.0 |     0.0 |
    | 03:00 |           1.0 |  0.5   0.5      0.0 |    1.0 |    0.0 |     0.5 |     0.5 |     0.0 |
    | 04:00 |           2.0 |  1.0   1.0      0.0 |    2.0 |    0.0 |     1.0 |     1.0 |     0.0 |
    | 05:00 |           3.0 |  1.5   0.5      1.0 |    3.0 |    0.0 |     1.5 |     0.5 |     1.0 |
    | 06:00 |           4.0 |  2.0   0.0      2.0 |    4.0 |    0.0 |     2.0 |     0.0 |     2.0 |
    | 07:00 |           5.0 |  2.5   0.0      4.0 |    5.0 |    0.0 |     2.5 |     0.0 |     4.0 |
    | 08:00 |           6.0 |  3.0   0.0      6.0 |    6.0 |    0.0 |     3.0 |     0.0 |     6.0 |
    | 09:00 |           7.0 |  3.5   0.0      8.0 |    7.0 |    0.0 |     3.5 |     0.0 |     8.0 |
    | 10:00 |           8.0 |  4.0   0.0     10.0 |    8.0 |    0.0 |     4.0 |     0.0 |    10.0 |

Adding discharge and keeping the water balance automatically works:

>>> rules(
...     PPolys(
...         output1=1,
...         output2=PPoly.from_data(xs=[0.0, 6.0, 7.0], ys=[0.0, 0.0, 1.0]),
...         output3=PPoly.from_data(xs=[-1.0, 0.0, 4.0, 6.0], ys=[-1.0, -1.0, 0.0, 0.0]),
...     )
... )

    >>> test("exch_branch_larsim_ex2")
    |  date | originalinput |            outputs | input1 | input2 | output1 | output2 | output3 |
    ----------------------------------------------------------------------------------------------
    | 00:00 |          -1.0 |  0.0  0.0     -1.0 |   -2.0 |    1.0 |     0.0 |     0.0 |    -1.0 |
    | 01:00 |          -1.0 |  0.0  0.0     -1.0 |   -1.0 |    0.0 |     0.0 |     0.0 |    -1.0 |
    | 02:00 |           0.0 |  1.0  0.0     -1.0 |    0.0 |    0.0 |     1.0 |     0.0 |    -1.0 |
    | 03:00 |           1.0 | 1.75  0.0    -0.75 |    1.0 |    0.0 |    1.75 |     0.0 |   -0.75 |
    | 04:00 |           2.0 |  2.5  0.0     -0.5 |    2.0 |    0.0 |     2.5 |     0.0 |    -0.5 |
    | 05:00 |           3.0 | 3.25  0.0    -0.25 |    3.0 |    0.0 |    3.25 |     0.0 |   -0.25 |
    | 06:00 |           4.0 |  4.0  0.0      0.0 |    4.0 |    0.0 |     4.0 |     0.0 |     0.0 |
    | 07:00 |           5.0 |  5.0  0.0      0.0 |    5.0 |    0.0 |     5.0 |     0.0 |     0.0 |
    | 08:00 |           6.0 |  6.0  0.0      0.0 |    6.0 |    0.0 |     6.0 |     0.0 |     0.0 |
    | 09:00 |           7.0 |  6.0  1.0      0.0 |    7.0 |    0.0 |     6.0 |     1.0 |     0.0 |
    | 10:00 |           8.0 |  6.0  2.0      0.0 |    8.0 |    0.0 |     6.0 |     2.0 |     0.0 |

Defining different rules for different seasons works:

>>> rule1 = PPolys(
...     output1=1,
...     output2=PPoly.from_data(xs=[-1.0, 0.0, 1.0], ys=[0.0, 0.0, 0.5]),
...     output3=PPoly.from_data(xs=[0.0, 1.0], ys=[0.0, 0.0]),
... )
>>> rule2 = PPolys(
...     output1=1,
...     output2=PPoly.from_data(xs=[0.0, 1.0], ys=[0.0, 0.0]),
...     output3=PPoly.from_data(xs=[-1.0, 0.0, 1.0], ys=[0.0, 0.0, 0.5]),
... )
>>> rules(
...     toy_01_01_00_30=rule1,
...     toy_01_01_04_30=rule1,
...     toy_01_01_09_30=rule2,
...     toy_12_31_23_30=rule2,
... )

At least I didn't find any inconsistencies in the show result tables.

Before tidying up, I will add the option to extract (or add) values from external time series instead of applying the PPolys rules.

We had a very thorough offline discussion. I hope the following list of the decisions we made is complete:

We will name the model developed so far HydPy-Exch-Branch-Rules, because it is only based on interpolation rules.
We will not add the option of specifying external exchange time series to it, because HydPy currently offers no way to store multiple input time series of the same sequence type for the same element instance in the same NetCDF, ASCII, or NumPy file.
We will implement the model HydPy-Exch-Branch-Ext instead, which expects a single exchange time series and falls back to a single interpolation rule if the time series is empty.
HydPy-Exch-Branch-Ext will generally work with the PPolys.REST idea (behind the scenes), so configuring it will be very easy and safe.
HydPy-Exch-Branch-Ext will work very much like the corresponding LARSIM features, but for now we ignore the option to define rules for specific time periods (meaning we will focus on seasonal patterns only).
HydPy-Exch-Branch-Rules will replace HydPy-Exch-Branch-HBV96 so that we will deprecate the latter.
For full compatibility, we must also add the HydPy-Exch-Branch-HBV96 parameters Delta and Minimum to HydPy-Exch-Branch-Rules. (We thought about renaming them to DeltaInput and MinimumInput for clarity, but - if I remember correctly - this was not a final decision.)
We will search for a way to simplify transferring HydPy-Exch-Branch-HBV96 configuration files to HydPy-Exch-Branch-Rules configuration files.
We will likely develop a separate submodel for the optional extraction of water from lakes and dams.

@parvumflumen

For VERZ and ERWV, LARSIM's documentation does not clarify what happens if the external time series (QE) requests more discharge than is currently flowing through the river (QR). I would expect one of the following behaviours:

QE stays as is, and QR is reduced to zero.
QE is reduced to QR, and QR is set to zero.
QE stays as is, and QR becomes negative.

After reading the available information on the same matter for VSPE, I think LARSIM favours option one because it considers the external time series as measured and certain. (@parvumflumen: Maybe you know better?)

However, violating the water balance equation is not always desired. Hence, we could define an option for choosing the preferred behaviour:

adjustextraction(ADJUST_OUTPUT)  # option 1
adjustextraction(ADJUST_EXCHANGE)  # option 2
adjustextraction(ADJUST_NOTHING)  # option 3

This parameter could also control the functionality of the fallback interpolation rule (if the given piecewise polynomials lead to such inconsistencies).

We considered naming the input sequence for providing the external time series Exchange. However, this name is already taken by the outlet sequence Exchange. My new suggestions:

Withdrawal
RequiredWithdrawal
Diversion
ExpectedDiversion

We considered naming the input sequence for providing the external time series Exchange. However, this name is already taken by the outlet sequence Exchange. My new suggestions:
* Withdrawal

* RequiredWithdrawal

* Diversion

* ExpectedDiversion

What about RequestedTransfer?, the other proposed names somehow carry a direction.

What about RequestedTransfer?, the other proposed names somehow carry a direction.

This was intentional because it clarifies the meaning of positive and negative values. But admittedly, negative withdrawal and negative diversion are a little unintuitive. RequestedTransfer is much better.

I also made some adjustments to the parameter names:

>>> targets(main="river", branch="diversion")
>>> fixwaterbalance(ADJUST_NOTHING)
>>> fallbackrules(PPoly.from_data(xs=[0.0, 2.0, 4.0, 6.0], ys=[0.0, 0.0, 2.0, 2.0]))

I think FixWaterBalance is much clearer than AdjustExtraction. It's current option names are:

ADJUST_NOTHING
ADJUST_BRANCH
ADJUST_RIVER

Instead of ADJUST_BRANCH, we could use ADJUST_TRANSFER, which is much clearer, especially after naming the input sequence RequestedTransferm. However, then ADJUST_RIVER doesn't sound right anymore. Maybe ADJUST_DISCHARGE?.

The name FallbackRules should clarify that the given rules only serve as a fallback. I first thought about reusing Rules, but a separate parameter allows simpler configuration and more thorough automatic validation.

Targets is a special parameter, that always has the integer value two. I needed it to clarify which outlet node serves which purpose. Alternatively, we could remove this parameter and use two different outlet sequence types, so this network-configuration aspect would be handled by Node.variable. This would be more consistent with how most models connect their sequences to nodes. However, models that branch their outflow to an arbitrary number of nodes definitely need such a model-internal selection mechanism. Hence, the proposed design is somewhat consistent with HydPy-Exch-Branch-Rules (and HydPy-Exch-Branch-HBV96).

Based on these considerations, my slightly modified suggestion:

>>> targets(stream="river", transfer="diversion")
>>> fixwaterbalance(ADJUST_NOTHING)  # ADJUST_TRANSFER / ADJUST_STREAM
>>> fallbackrules(PPoly.from_data(xs=[0.0, 2.0, 4.0, 6.0], ys=[0.0, 0.0, 2.0, 2.0]))

Deprecate `exch_branch_hbv96` while adding the transition method `con…

Add the branching model `exch_branch_io` (see issue #201).

Add the branching model `exch_branch_rules` (see issue #201).

Deprecate `exch_branch_hbv96` while adding the transition method `con…

Add the branching model `exch_branch_io` (see issue #201).

Add the branching model `exch_branch_rules` (see issue #201).

Deprecate `exch_branch_hbv96` while adding the transition method `con…

Add the branching model `exch_branch_io` (see issue #201).

Add the branching model `exch_branch_rules` (see issue #201).

Deprecate `exch_branch_hbv96` while adding the transition method `con…

The changes discussed here are now included in the master branch and in the documentation. The following short examples should cover our final (naming) decisions:

>>> from hydpy.models.exch_branch_io import *
>>> parameterstep()
>>> targets(stream="river", transfer="diversion")
>>> inputs.requestedtransfer.series = 6 * [2.0] + 5 * [-2.0]
>>> minstream(1.0)
>>> maxstream(8.0)
>>> keepwaterbalance(True)
>>> flowtransferrules(
...     toy_01_01_12=PPoly(xs=[0.0, 2.0, 4.0, 6.0], ys=[0.0, 0.0, 1.0, 1.0]),
...     toy_01_09_12=PPoly(xs=[0.0, 2.0, 4.0, 6.0], ys=[0.0, 0.0, 1.0, 1.0]),
...     toy_01_11_12=PPoly(xs=[0.0, 2.0, 4.0, 6.0], ys=[0.0, 0.0, 2.0, 2.0]),
... )

>>> from hydpy.models.exch_branch_rules import *
>>> parameterstep()
>>> rules(
...     PPolys(
...         output1=PPolys.REST,
...         output2=PPoly(xs=[0.0, 6.0, 7.0], ys=[0.0, 0.0, 1.0]),
...         output3=PPoly(xs=[-1.0, 0.0, 4.0, 6.0], ys=[-1.0, -1.0, 0.0, 0.0]),
...     )
... )

exch_branch_hbv96 is now deprecated and can be automatically converted to exch_branch_rules with method convert_2_exch_branch_rules. We nevertheless did not add the parameters Delta and Minimum to exch_branch_rules, because they can be simply emulated by piecewise interpolation rules (as shown in the method convert_2_exch_branch_rules and the integration test input-modifying rules).

We switched from the name FallbackRules to FlowTransferRules because exch_water_transfer_io requires a similar but not identical fallback parameter. Since we cannot use the same name twice, we chose FlowTransferRules and VolTransferRules, which clarify the different input types.

Add the branching model `exch_branch_io` (see issue #201).

Add the branching model `exch_branch_rules` (see issue #201).

Deprecate `exch_branch_hbv96` while adding the transition method `con…

Navigation Menu

Repository navigation

Support more of LARSIM's branching features #201

`Rules` as a "special" parameter

`Rules` subclassed from SeasonalInterpolator

`Rules` subclassed from SeasonalInterpolator

10 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Support more of LARSIM's branching features #201

Description

Activity

tyralla commented on Sep 21, 2026

tyralla commented on Sep 22, 2026

tyralla commented on Sep 22, 2026

Rules as a "special" parameter

Rules subclassed from SeasonalInterpolator

tyralla commented on Sep 24, 2026

Rules subclassed from SeasonalInterpolator

tyralla commented on Sep 24, 2026

tyralla commented on Sep 25, 2026

tyralla commented on Sep 25, 2026

tyralla commented on Sep 26, 2026

tyralla commented on Sep 26, 2026

BGWKlein commented on Sep 28, 2026

tyralla commented on Sep 28, 2026

10 remaining items

tyralla commented on Oct 11, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

`Rules` as a "special" parameter

`Rules` subclassed from SeasonalInterpolator

`Rules` subclassed from SeasonalInterpolator