Author: Marco Dal Molin
Collaborators: Fabrizio Fenicia, Dmitri Kavetski
This example will cover the following aspects:
By clicking on the items you will be redirected to the documentation page (when available) that explains the arguments in a more detailed way.
The following aspects are already covered in other examples:
For this reason we will put only the code need, without further explanations. You can check the other examples by clicking on the items above.
Examples of SuperflexPy are created and distributed using Jupyter notebooks because they enable to conjugate runnable code with text explanations. We have decided to not provide the content of the notebooks as script because this would mean duplicating the same content in multiple places, bringing to maintainability problems.
If the user wants to download the content of this notebook as a python script, it is possible following the steps:
The result is a Python script that contains all the markdown text as comment and the code contained inside the cells.
The execution of the code in this example relies on the following external libraries:
We assume that those libraries are already installed together with the latest version of SuperflexPy. Keep in mind that not all the libraries listed above are strictly needed to execute SuperflexPy and, therefore, they are not listed as requirements of SuperflexPy.
import numpy as np
from datetime import datetime
In this example we will use the Unit as maximum complexity level of the model. The same operations illustrated in this example can be done similarly with any SuperflexPy component (i.e. Element, Unit, Node, Network).
The model is imported and initialized in the following cell.
from superflexpy.implementation.root_finders.pegasus import PegasusPython
from superflexpy.implementation.numerical_approximators.implicit_euler import ImplicitEulerPython
from superflexpy.implementation.elements.hbv import PowerReservoir
from superflexpy.framework.unit import Unit
# Initialize numercal routines
root_finder = PegasusPython()
numeric_approximator = ImplicitEulerPython(root_finder=root_finder)
# Initialize the elements
pr_1 = PowerReservoir(
parameters={'k': 0.1, 'alpha': 1.0},
states={'S0': 10.0},
approximation=numeric_approximator,
id='PR-1'
)
pr_2 = PowerReservoir(
parameters={'k': 0.01, 'alpha': 2.0},
states={'S0': 1.0},
approximation=numeric_approximator,
id='PR-2'
)
# Initialize the Unit
model = Unit(
layers=[
[pr_1],
[pr_2]
],
id='model'
)
This involves:
This is done in the following cell.
# Fix the seed
SEED = 2
rng = np.random.RandomState(seed=SEED)
# Create the precipitation input
P = np.zeros(100)
P[:10] = rng.randint(10, size=10)
P[25:30] = rng.randint(20, size=5)
P[40:60] = rng.randint(5, size=20)
P[80:83] = rng.randint(30, 50, size=3)
# Assign the input
model.set_input([P])
# Set the time step
model.set_timestep(1.0)
This example will, actually, not show how to calibrate the model to data but it will illustrate how to perform all the steps that are usually needed to calibrate a model. The steps are:
Steps 1 to 4 are usually repeated several times until reaching satisfactory values for the parameters.
We want to calibrate only 2 out of 4 parameters: model_PR-1_k
and model_PR-2_k
. Normally, the calibration library will provide the value of these parameters. In this case, we will generate them randomly since the goal is to show the process and not to actually calibrate the model.
# Generate the parameters randomly
parameters = rng.rand(2)
# Construct the dictionary for the set_parameters method
parameters_to_set = {
'model_PR-1_k': parameters[0],
'model_PR-2_k': parameters[1],
}
This can be done using the set_parameters
method
model.set_parameters(parameters=parameters_to_set)
To avoid errors related to previous runs of the model, it is good practice to reset the states to their initialization value.
model.reset_states()
No the model can be run using the get_output
method.
output = model.get_output()
During calibration (or uncertainty quantification) the steps illustrated above can be repeated several hundreds (if not thousands) of times. For this reason, using the pure Python implementation of the numerical solvers may not be the most efficient thing to do.
We will here show the advantage of using the numba implementation. Further details can be found here.
This function will performs the steps seen in "STEP 3".
def calibrate(model, rng, num_iterations):
for _ in range(num_iterations):
parameters = rng.rand(2)
parameters_to_set = {
'model_PR-1_k': parameters[0],
'model_PR-2_k': parameters[1],
}
model.set_parameters(parameters=parameters_to_set)
model.reset_states()
output = model.get_output()
The only thing that need to be change in order to use the model with numba is the choice of the root finder and of the numerical approximator to use.
from superflexpy.implementation.root_finders.pegasus import PegasusNumba
from superflexpy.implementation.numerical_approximators.implicit_euler import ImplicitEulerNumba
# Initialize numercal routines
root_finder_numba = PegasusNumba()
numeric_approximator_numba = ImplicitEulerNumba(root_finder=root_finder_numba)
# Initialize the elements
pr_1_numba = PowerReservoir(
parameters={'k': 0.1, 'alpha': 1.0},
states={'S0': 10.0},
approximation=numeric_approximator_numba,
id='PR-1'
)
pr_2_numba = PowerReservoir(
parameters={'k': 0.01, 'alpha': 2.0},
states={'S0': 1.0},
approximation=numeric_approximator_numba,
id='PR-2'
)
# Initialize the Unit
model_numba = Unit(
layers=[
[pr_1_numba],
[pr_2_numba]
],
id='model'
)
# Assign the input
model_numba.set_input([P])
# Set the time step
model_numba.set_timestep(1.0)
We will simulate a calibration that requires 1000 model runs to calibrate.
start = datetime.now()
calibrate(model=model, rng=rng, num_iterations=1000)
end = datetime.now()
print('Calibration time: {}'.format(end - start))
Calibration time: 0:00:27.209784
We will simulate a calibration that requires 1000 model runs to calibrate.
start = datetime.now()
calibrate(model=model_numba, rng=rng, num_iterations=1000)
end = datetime.now()
print('Calibration time: {}'.format(end - start))
Calibration time: 0:00:05.353659
Although running times depends on the hardware used, numba models can be up to 30 times faster than pure Python solutions.
Note that, when using numba for the first time, most of the time is used to compile the methods to machine code (in the machine used when writing the example, compilation takes roughly 1 second while running the model once is in the order of $10^{-3}$ seconds). Compilation needs to be done only once as long as data types do not change and, therefore, it is more efficient to use numba when the model needs to be run several times, like for calibration.