Multiple run#
Here, we aim to compute multiple forward runs of a single Model object (in parallel) in order to retrieve the cost (and discharge values) for each parameters set tested. This method is often used to perform global sensitivity analysis.
To get started, open a Python interface:
python3
Imports#
In [1]: import smash
In [2]: import numpy as np
In [3]: import time
Model object creation#
First, you need to create a smash.Model object.
For this case, we will use the Lez dataset.
Load the setup and mesh dictionaries using the smash.load_dataset() method and create the smash.Model object.
In [4]: setup, mesh = smash.load_dataset("Lez")
In [5]: model = smash.Model(setup, mesh)
Generation of the sample#
In order to compute several forward runs of the smash.Model object, we first need to generate a sample of parameters.
We define a custom sample with random sets of parameters using the smash.generate_samples() method.
First, we refer to the smash.Model.get_bound_constraints() method to obtain some information about the Model parameters:
In [6]: model.get_bound_constraints()
Out[6]:
{'num_vars': 4,
'names': ['cp', 'cft', 'exc', 'lr'],
'bounds': [[9.999999974752427e-07, 1000.0],
[9.999999974752427e-07, 1000.0],
[-50.0, 50.0],
[9.999999974752427e-07, 1000.0]]}
Then, we will define a problem based on the previous information by adjusting the bounds of the parameters:
In [7]: problem = {
...: "num_vars": 4,
...: "names": ["cp", "cft", "exc", "lr"],
...: "bounds": [[1, 2000], [1, 1000], [-20, 5], [1, 1000]]
...: }
...:
Once the problem is defined we can generate the sample. Here we define a sample of 400 random sets of Model parameters:
In [8]: sample = smash.generate_samples(problem, n=400, random_state=1)
Note
random_state argument is used to set the random seed.
Visualization of the sample#
sample is a SampleResult object with several methods including the SampleResult.slice() method.
This method allows to create slices of the original sample (sub-sample). We will use this method to quickly visualize the sample.
Visualization of the 10 first sets:
In [9]: sample.slice(10)
Out[9]:
cft: array([959.47488674, 804.15693061, 33.29074356, 709.67786398,
465.53648008, 947.60139245, 222.21130159, 267.80494766,
82.39249091, 429.19020984])
cp: array([8.34626987e+02, 1.44092866e+03, 1.22863526e+00, 6.05362813e+02,
2.94365026e+02, 1.85584851e+02, 3.73334163e+02, 6.91775893e+02,
7.94138181e+02, 1.07809465e+03])
exc: array([ 2.20966086, -15.39038894, -5.36628853, 2.45512615,
-8.84706952, 3.04670764, -13.02522774, -4.77922065,
-2.93865726, -14.29485667])
generator: 'uniform'
lr: array([822.02752539, 628.61135753, 956.34520825, 590.33320395,
198.63495186, 429.90139189, 337.4092527 , 991.9530669 ,
380.85625962, 992.70817279])
n_sample: 10
We can also visualize sets between start and end index:
In [10]: sample.slice(start=10, end=20)
Out[10]:
cft: array([109.90974526, 634.15297385, 803.16027414, 697.10369562,
766.44516969, 343.11166572, 846.00563154, 429.34000622,
824.18586054, 626.86966247])
cp: array([ 838.96983429, 1370.75378129, 409.70004721, 1756.35675535,
55.7477988 , 1341.26455285, 835.19229993, 1117.82096706,
281.63349025, 397.00487668])
exc: array([-19.65581219, -9.58190094, 3.4620473 , -11.42429725,
-0.50639261, -15.63159218, -11.45117909, -16.38505698,
-2.08072965, -2.51730941])
generator: 'uniform'
lr: array([519.28544648, 173.10742178, 75.49949718, 370.93538808,
124.29544095, 634.88419758, 414.50592588, 990.92471689,
930.02566803, 150.08679628])
n_sample: 10
Conversion of the sample#
Two other methods can be used to convert this object to pandas.DataFrame or numpy.ndarray. The SampleResult.to_dataframe() method and
SampleResult.to_numpy()
In [11]: sample.to_dataframe()
Out[11]:
cp cft exc lr
0 834.626987 959.474887 2.209661 822.027525
1 1440.928662 804.156931 -15.390389 628.611358
2 1.228635 33.290744 -5.366289 956.345208
3 605.362813 709.677864 2.455126 590.333204
4 294.365026 465.536480 -8.847070 198.634952
.. ... ... ... ...
395 1687.837017 808.469736 -3.322030 595.252286
396 762.651232 295.993506 -11.391520 878.782315
397 1499.966763 544.577262 -16.215611 255.346155
398 1022.771815 488.433570 -4.115726 766.422177
399 1082.362658 855.501056 1.194701 941.072254
[400 rows x 4 columns]
# axis=-1 to stack along columns
In [12]: sample.to_numpy(axis=-1)
Out[12]:
array([[ 8.34626987e+02, 9.59474887e+02, 2.20966086e+00,
8.22027525e+02],
[ 1.44092866e+03, 8.04156931e+02, -1.53903889e+01,
6.28611358e+02],
[ 1.22863526e+00, 3.32907436e+01, -5.36628853e+00,
9.56345208e+02],
...,
[ 1.49996676e+03, 5.44577262e+02, -1.62156108e+01,
2.55346155e+02],
[ 1.02277182e+03, 4.88433570e+02, -4.11572596e+00,
7.66422177e+02],
[ 1.08236266e+03, 8.55501056e+02, 1.19470109e+00,
9.41072254e+02]])
Computation of the multiple runs#
Once the sample is generated, the mutliple runs can be simulated using the Model.multiple_run method.
Multiple runs with default setup#
Here we will compute the multiple runs with the default setup. That is, returning the value of the cost function for each set using 1 CPU.
The default cost function is a nse calculated on the most downstream gauge.
In [13]: mtprr = model.multiple_run(sample)
</> Multiple Run Model
Nsample: 400
Ncpu: 1
Jobs function: [ nse ]
wJobs function: [ 1.0 ]
Np: 4 [ cp cft exc lr ]
Ns: 0 [ ]
Ng: 1 [ Y3204040 ]
wg: 1 [ 1.0 ]
mtprr is a MultipleRunResult object which has a cost attribute. We can visualize the 10 first cost values returned.
In [14]: mtprr.cost[0:10]
Out[14]:
array([1.349332 , 1.3831991, 0.5622262, 1.2688785, 0.8765431, 0.9696094,
0.820316 , 1.2330167, 0.9808377, 1.372884 ], dtype=float32)
Multiple runs with simulated discharge returned#
We can return the simulated discharge for each set by setting True to the return_qsim argument.
In [15]: mtprr_qsim = model.multiple_run(sample, return_qsim=True)
</> Multiple Run Model
Nsample: 400
Ncpu: 1
Jobs function: [ nse ]
wJobs function: [ 1.0 ]
Np: 4 [ cp cft exc lr ]
Ns: 0 [ ]
Ng: 1 [ Y3204040 ]
wg: 1 [ 1.0 ]
This time, mtprr_qsim is a MultipleRunResult object which has a cost and qsim attribute. We can visualize the 10 first cost and simulated
discharge (on the most downstream gauge) values returned.
In [16]: mtprr_qsim.cost[0:10]
Out[16]:
array([1.349332 , 1.3831991, 0.5622262, 1.2688785, 0.8765431, 0.9696094,
0.820316 , 1.2330167, 0.9808377, 1.372884 ], dtype=float32)
# The shape correspond to (number of gauges, number of time steps, number of sets)
In [17]: mtprr_qsim.qsim.shape
Out[17]: (3, 364, 400)
In [18]: mtprr_qsim.qsim[0,:,0:10]
Out[18]:
array([[8.3887035e-06, 1.3336102e-05, 6.4033538e-06, ..., 6.0160660e-06,
2.5853311e-05, 6.0082957e-06],
[7.0552037e-06, 8.8686575e-06, 5.6740519e-06, ..., 5.3902600e-06,
5.5817004e-06, 5.3844692e-06],
[5.6783142e-06, 4.8928027e-06, 4.8806255e-06, ..., 4.7030931e-06,
9.4816346e-07, 4.6993473e-06],
...,
[7.3203817e-02, 1.4144397e-02, 8.5565358e-02, ..., 3.1998256e-01,
1.3792305e-01, 1.1160661e-01],
[7.6349042e-02, 1.4939376e-02, 8.5804544e-02, ..., 3.1521481e-01,
1.3526769e-01, 1.1082144e-01],
[7.8357331e-02, 1.2541366e-02, 8.8809021e-02, ..., 3.1069714e-01,
1.3218015e-01, 1.1008787e-01]], dtype=float32)
Multiple runs with customize cost options#
We can also change the cost option arguments such as the kind of objective function (jobs_fun argument) and on which gauge we want to compute the cost (gauge argument).
In [19]: mtprr_cost_opt = model.multiple_run(sample, jobs_fun="kge", gauge="all")
</> Multiple Run Model
Nsample: 400
Ncpu: 1
Jobs function: [ kge ]
wJobs function: [ 1.0 ]
Np: 4 [ cp cft exc lr ]
Ns: 0 [ ]
Ng: 3 [ Y3204040 Y3204030 Y3204010 ]
wg: 3 [ 0.33333334 0.33333334 0.33333334 ]
We can visualize the 10 first cost values returned.
In [20]: mtprr_cost_opt.cost[0:10]
Out[20]:
array([1.3900341 , 1.4084475 , 0.39728892, 1.3754207 , 1.0251681 ,
1.1090906 , 1.0353888 , 1.3561141 , 1.1567951 , 1.4364334 ],
dtype=float32)
Multiple runs in parallel#
At the moment, all the previous runs were done sequentially. we can save computation time by parallelizing the runs on several CPUs by assigning a value to the ncpu argument.
We will used the time library, previously imported, to retrieve the computation time and compare it between a sequential and parallel run.
In [21]: start_seq = time.time()
In [22]: mtprr_seq = model.multiple_run(sample)
</> Multiple Run Model
Nsample: 400
Ncpu: 1
Jobs function: [ nse ]
wJobs function: [ 1.0 ]
Np: 4 [ cp cft exc lr ]
Ns: 0 [ ]
Ng: 1 [ Y3204040 ]
wg: 1 [ 1.0 ]
In [23]: time_seq = time.time() - start_seq
In [24]: start_prl = time.time()
In [25]: mtprr_prl = model.multiple_run(sample, ncpu=4)
</> Multiple Run Model
Nsample: 400
Ncpu: 4
Jobs function: [ nse ]
wJobs function: [ 1.0 ]
Np: 4 [ cp cft exc lr ]
Ns: 0 [ ]
Ng: 1 [ Y3204040 ]
wg: 1 [ 1.0 ]
In [26]: time_prl = time.time() - start_prl
We can now visualize that we have the same data between a sequential and parallel run and the difference of time computation in seconds.
In [27]: mtprr_seq.cost[0:10]
Out[27]:
array([1.349332 , 1.3831991, 0.5622262, 1.2688785, 0.8765431, 0.9696094,
0.820316 , 1.2330167, 0.9808377, 1.372884 ], dtype=float32)
In [28]: mtprr_prl.cost[0:10]
Out[28]:
array([1.349332 , 1.3831991, 0.5622262, 1.2688785, 0.8765431, 0.9696094,
0.820316 , 1.2330167, 0.9808377, 1.372884 ], dtype=float32)
In [29]: time_seq, time_prl
Out[29]: (3.5197415351867676, 1.1425659656524658)
Multiple runs with iterator#
In case you want to retrieve the simulated discharge for each set of the sample, it may be impossible to store the entire data in memory.
As a reminder, the simulated discharge data is stored as single precision floating point (float32) in an array of size \(ng \times nt \times ns\) where \(ng\) is the number of gauges,
\(nt\) the number of time steps and \(ns\) the number of sets in the sample.
In our case, it results in an array of size \(n=3 \times 364 \times 400 = 436 800\). One way to get around this problem is to iterate on sub-samples and to serialize each output to
free the memory.
To perform this we will use the SampleResult.iterslice() method to iterate on the sample.
This method as a by argument, set to 1 by default, which allows the user to choose the size of the slice (sub-sample) to iterate on.
We will fix the by argument to 100 to perform 4 iterations on a sample of size 400.
In this way, we will divide by 4 the memory space taken by the array of simulated discharge.
In [30]: for i, slc in enumerate(sample.iterslice(by=100)):
....: mtprr = model.multiple_run(slc, ncpu=4, return_qsim=True)
....: smash.save_model_ddt(
....: model,
....: f"res_mtprr_slc_{i}.hdf5",
....: sub_data={"mtprr_cost": mtprr.cost, "mtprr_qsim": mtprr.qsim}
....: )
....:
Note
We used the smash.save_model_ddt() method to serialize the result on each iteration by filling in the sub_data argument.