From bd8bb8407358a50876ddc41692dca816de844011 Mon Sep 17 00:00:00 2001 From: thouska Date: Fri, 25 Sep 2020 14:08:36 +0200 Subject: [PATCH] Introduction of new documentation slides #259 --- docs/Advanced_hints.md | 24 +-- ...ayesian_uncertainty_analysis_with_DREAM.md | 128 +++++++++++++ ...brating a hydrological model with DREAM.md | 174 ------------------ docs/Calibration_with_SCE-UA.md | 88 +++++++++ docs/How_to_link_a_model_to_SPOTPY.md | 126 +++++++++++++ docs/Sensitivity_analysis_with_FAST.md | 63 +++++++ docs/index.md | 6 +- mkdocs.yml | 5 +- spotpy/algorithms/sceua.py | 2 +- spotpy/analyser.py | 2 +- spotpy/examples/spot_setup_hymod_python.py | 2 +- spotpy/examples/tutorial_fast_hymod.py | 11 +- spotpy/examples/tutorial_sceua_hymod.py | 128 +++---------- 13 files changed, 441 insertions(+), 318 deletions(-) create mode 100644 docs/Bayesian_uncertainty_analysis_with_DREAM.md delete mode 100644 docs/Calibrating a hydrological model with DREAM.md create mode 100644 docs/Calibration_with_SCE-UA.md create mode 100644 docs/How_to_link_a_model_to_SPOTPY.md create mode 100644 docs/Sensitivity_analysis_with_FAST.md diff --git a/docs/Advanced_hints.md b/docs/Advanced_hints.md index da74b6c3..a86e1166 100644 --- a/docs/Advanced_hints.md +++ b/docs/Advanced_hints.md @@ -20,7 +20,7 @@ To tell SPOTPY to use MPI, just give this information to the sampler: sampler = spotpy.algorithms.sceua(spotpy_setup,() dbname='RosenSCEUA', dbformat='csv', parallel='mpi') sampler.sample(10000) - + Now save this file and start it from a console: `mpirun -c 20 your_script.py`, where 20 is the number of cores you want to use. This should give you and speed of neerly 20 times compared with the standard sequential sampling. @@ -38,29 +38,7 @@ DE-MCz will parallelize the selcted number of chains [terBrack and Vrugt, 2008]( Th algorithms `MLE`, `MCMC` and `SA` can not run in parallel. -## FAST - Sensitivity analysis -SPOTPY gives you the opportunity to start a sensitivity analysis of your model. In this case, we included a global sensitivity analysis called "Extended FAST" based on -Saltelli et al. (1999). This is besides the Sobol´ sensitivity test the only algorithm available that is taking parameter interaction into account. - -The algorithm will tell you, how sensitive your parameters are on whatever is given back by your objective function. Before you start to sample, you should know how how many -iterations you need to get an reliable information about your parameter. The number of iteration can be calculate after [Henkel et al. GLOBAL SENSITIVITY ANALYSIS OF NONLINEAR MATHEMATICAL MODELS - AN -IMPLEMENTATION OF TWO COMPLEMENTING VARIANCE-BASED ALGORITHMS, 2012] (https://www.informs-sim.org/wsc12papers/includes/files/con308.pdf): - -$$N = (1+4M^2(1+(k-2)d))k$$ - -with N = needed parameter iterations, M= inference factor (SPOTPY default M=4) and d= frequency step (SPOTPY default d=2) and k as the number of parameters of your model. -You can start the simulation with - - sampler = spotpy.algorithms.fast(spotpy_setup,() dbname='Fast_sensitivity', dbformat='csv') - -and you can analyse your results with - - results = sampler.get_data() - analyser.plot_fast_sensitivity(results,number_of_sensitiv_pars=5) - -This will show you a Plot with the total Sensitivity index of all your parameters and in this case the five most sensitive parameters (can be adjusted). - ## Plotting time If you want create plots out of your samples and you don't want to sample your results again do something like this: diff --git a/docs/Bayesian_uncertainty_analysis_with_DREAM.md b/docs/Bayesian_uncertainty_analysis_with_DREAM.md new file mode 100644 index 00000000..c0a4d104 --- /dev/null +++ b/docs/Bayesian_uncertainty_analysis_with_DREAM.md @@ -0,0 +1,128 @@ +# Bayesian uncertainty analysis of HYMOD with DREAM + +This chapter shows you, how to calibrate an external hydrological model (HYMOD) with SPOTPY. + +We use the previously created hydrological model HYMOD spotpy_setup class as an [example](https://github.com/thouska/spotpy/blob/master/spotpy/examples/spot_setup_hymod_python.py), +to perform a parameter uncertainty analysis and Bayesian calibration with the Differential Evolution Adaptive Metropolis (DREAM) algorithm. +For detailed information about the underlying theory, have a look at the [Vrugt (2016)](https://doi.org/10.1016/j.envsoft.2015.08.013 "Vrugt (2016)"). + + +First some relevant functions are imported Hymod example: + + import numpy as np + import spotpy + import matplotlib.pyplot as plt + from spotpy.likelihoods import gaussianLikelihoodMeasErrorOut as GausianLike + from spotpy.analyser import plot_parameter_trace + from spotpy.analyser import plot_posterior_parameter_histogram + +Further we need the spotpy_setup class, which links the model to spotpy: + + from spotpy.examples.spot_setup_hymod_python import spot_setup + +It is important that this function is initialized before it is further used by SPOTPY: + + spot_setup=spot_setup() + +In this special case, we want to change the implemented rmse objective function with a likelihood function, which is compatible with the Bayesian +calibration approach, used by the DREAM sampler: + + # Initialize the Hymod example (will only work on Windows systems) + spot_setup=spot_setup(GausianLike) + + +## Sample with DREAM + +Now we create a sampler, by using one of the implemented algorithms in SPOTPY. +The algorithm needs the inititalized spotpy setup class, and wants you to define a database name and database format. Here we create a DREAM_hymod.csv file, +which will contain the returned likelihood function, the coresponding parameter set, the model results and a chain ID (Dream is starting different independent +optimization chains, which need to be in line, to receive robust results): + + sampler=spotpy.algorithms.dream(spot_setup, dbname='DREAM_hymod', dbformat='csv') + +To actually start the algorithm, spotpy needs some further details, like the maximum allowed number of repetitions, the number of chains used by dream (default = 5) and set the Gelman-Rubin convergence limit (default 1.2). +We further allow 100 runs after convergence is achieved: + + #Select number of maximum repetitions + rep=5000 + + # Select five chains and set the Gelman-Rubin convergence limit + nChains = 4 + convergence_limit = 1.2 + + # Other possible settings to modify the DREAM algorithm, for details see Vrugt (2016) + nCr = 3 + eps = 10e-6 + runs_after_convergence = 100 + acceptance_test_option = 6 + +We start the sampler and collect the gained r_hat convergence values after the sampling: + + r_hat = sampler.sample(rep, nChains, nCr, eps, convergence_limit) + + +## Access the results +All gained results can be accessed from the SPOTPY csv-database: + + results = spotpy.analyser.load_csv_results('DREAM_hymod') + +These results are structured as a numpy array. Accordingly, you can access the different columns by using simple Python code, +e.g. to access all the simulations: + + fields=[word for word in results.dtype.names if word.startswith('sim')] + print(results[fields) + +## Plot model uncertainty +For the analysis we provide some examples how to plor the data. +If you want to see the remaining posterior model uncertainty: + + fig= plt.figure(figsize=(16,9)) + ax = plt.subplot(1,1,1) + q5,q25,q75,q95=[],[],[],[] + for field in fields: + q5.append(np.percentile(results[field][-100:-1],2.5)) + q95.append(np.percentile(results[field][-100:-1],97.5)) + ax.plot(q5,color='dimgrey',linestyle='solid') + ax.plot(q95,color='dimgrey',linestyle='solid') + ax.fill_between(np.arange(0,len(q5),1),list(q5),list(q95),facecolor='dimgrey',zorder=0, + linewidth=0,label='parameter uncertainty') + ax.plot(spot_setup.evaluation(),'r.',label='data') + ax.set_ylim(-50,450) + ax.set_xlim(0,729) + ax.legend() + fig.savefig('python_hymod.png',dpi=300) + +![Posterior model uncertainty](../img/python_hymod_simulation.png) +*Figure 1: Posterior model uncertainty of HYMOD.* + +## Plot convergence diagnostic +If you want to check the convergence of the DREAM algorithm: + + spotpy.analyser.plot_gelman_rubin(results, r_hat, fig_name='python_hymod_convergence.png') + +![Convergence diagnostic](../img/python_hymod_convergence.png) + +*Figure 2: Gelman-Rubin onvergence diagnostic of DREAM results.* + + +## Plot parameter uncertainty +Or if you want to check the posterior parameter distribution: + + parameters = spotpy.parameter.get_parameters_array(spot_setup) + + fig, ax = plt.subplots(nrows=5, ncols=2) + for par_id in range(len(parameters)): + plot_parameter_trace(ax[par_id][0], results, parameters[par_id]) + plot_posterior_parameter_histogram(ax[par_id][1], results, parameters[par_id]) + + ax[-1][0].set_xlabel('Iterations') + ax[-1][1].set_xlabel('Parameter range') + + plt.show() + fig.savefig('hymod_parameters.png',dpi=300) + +![Parameter distribution](../img/python_hymod_parameters.png) + +*Figure 3: Posterior parameter distribution of HYMOD.* + +The corresponding code is available for download [here](https://github.com/thouska/spotpy/blob/master/spotpy/examples/tutorial_dream_hymod.py). \ No newline at end of file diff --git a/docs/Calibrating a hydrological model with DREAM.md b/docs/Calibrating a hydrological model with DREAM.md deleted file mode 100644 index d85f79e7..00000000 --- a/docs/Calibrating a hydrological model with DREAM.md +++ /dev/null @@ -1,174 +0,0 @@ -# Calibration of HYMOD with DREAM - -This chapter shows you, how to link the external hydrological model HYMOD with SPOTPY (works only on Windows systems). - -We use the hydrological model HYMOD as an example, to calibrate it with the Differential Evolution Adaptive Metropolis (DREAM) algorithm. -For detailed information about the underlying theorie, have a look at the [Vrugt (2016)](https://doi.org/10.1016/j.envsoft.2015.08.013 "Vrugt (2016)"). -The SPOTPY package comes with an example which is desgined to help you to set up your own research project. - -First we need to import some functions we will need in the following: - - from spotpy.parameter import Uniform - from spotpy.examples.hymod_python.hymod import hymod - import os - -Now, we can setup the model within a spot setup class. The model needs some meteorological input data and five parameters to estimate discharge: - -## Connect HYMOD with SPOTPY - -Here we use class parameter, to initialize the parameter for our model. - - class spot_setup(object): - cmax = Uniform(low=1.0 , high=500, optguess=412.33) - bexp = Uniform(low=0.1 , high=2.0, optguess=0.1725) - alpha = Uniform(low=0.1 , high=0.99, optguess=0.8127) - Ks = Uniform(low=0.001 , high=0.10, optguess=0.0404) - Kq = Uniform(low=0.1 , high=0.99, optguess=0.5592) - - def __init__(self, obj_func=None): - #Just a way to keep this example flexible and applicable to various examples - self.obj_func = obj_func - #Transform [mm/day] into [l s-1], where 1.783 is the catchment area - self.Factor = 1.783 * 1000 * 1000 / (60 * 60 * 24) - self.PET,self.Precip = [], [] - self.date,self.trueObs = [], [] - #Find Path to Hymod on users system - self.owd = os.path.dirname(os.path.realpath(__file__)) - self.hymod_path = self.owd+os.sep+'hymod_python' - #Load Observation data from file - climatefile = open(self.hymod_path+os.sep+'hymod_input.csv', 'r') - headerline = climatefile.readline()[:-1] - - #Read model forcing in working storage (this is done only ones) - if ';' in headerline: - self.delimiter = ';' - else: - self.delimiter = ',' - self.header = headerline.split(self.delimiter) - for line in climatefile: - values = line.strip().split(self.delimiter) - self.date.append(str(values[0])) - self.Precip.append(float(values[1])) - self.PET.append(float(values[2])) - self.trueObs.append(float(values[3])) - climatefile.close() - -We use the simulation function to write one random parameter set into a parameter file, like it is needed for the HYMOD model, -start the model and read the model discharge output data: - - def simulation(self,x): - #Here the model is actualy startet with one paramter combination that it gets from spotpy for each time the model is called - data = hymod(self.Precip, self.PET, x[0], x[1], x[2], x[3], x[4]) - sim=[] - for val in data: - sim.append(val*self.Factor) - #The first year of simulation data is ignored (warm-up) - return sim[366:] - -And in a last step, we compare the observed and the simulated data. Here we set up the setup class in a way, that it can receive any -likein the SPOTPY package. Please mind that the selection of the Likelihood highly influences the results gained with this algorithm: - - def evaluation(self): - #The first year of simulation data is ignored (warm-up) - return self.trueObs[366:] - - def objectivefunction(self,simulation,evaluation, params=None): - like = spotpy.likelihoods.gaussianLikelihoodMeasErrorOut(evaluation,simulation) - return like - -## Sample with DREAM - -Now we can initialize the Hymod example: - - spot_setup=spot_setup() - -Create the Dream sampler of spotpy, alt_objfun is set to None to force SPOTPY -to jump into the def objectivefunction in the spot_setup class (default is -spotpy.objectivefunctions.log_p). Results are saved in a DREAM_hymod.csv file: - - sampler=spotpy.algorithms.dream(spot_setup, dbname='DREAM_hymod', dbformat='csv') - -Select number of maximum repetitions, the number of chains used by dream (default = 5) and set the Gelman-Rubin convergence limit (default 1.2). -We further allow 100 runs after convergence is achieved: - - #Select number of maximum repetitions - rep=5000 - - # Select five chains and set the Gelman-Rubin convergence limit - nChains = 4 - convergence_limit = 1.2 - - # Other possible settings to modify the DREAM algorithm, for details see Vrugt (2016) - nCr = 3 - eps = 10e-6 - runs_after_convergence = 100 - acceptance_test_option = 6 - -We start the sampler and collect the gained r_hat convergence values after the sampling: - - r_hat = sampler.sample(rep, nChains, nCr, eps, convergence_limit) - - -## Access the results -All other results can be accessed from the SPOTPY csv-database: - - results = spotpy.analyser.load_csv_results('DREAM_hymod') - -These results are structured as a numpy array. Accordingly, you can access the different columns by using simple Python code, -e.g. to access all the simulations: - - fields=[word for word in results.dtype.names if word.startswith('sim')] - print(results[fields) - -## Plot model uncertainty -For the analysis we provide some examples how to plor the data. -If you want to see the remaining posterior model uncertainty: - - fig= plt.figure(figsize=(16,9)) - ax = plt.subplot(1,1,1) - q5,q25,q75,q95=[],[],[],[] - for field in fields: - q5.append(np.percentile(results[field][-100:-1],2.5)) - q95.append(np.percentile(results[field][-100:-1],97.5)) - ax.plot(q5,color='dimgrey',linestyle='solid') - ax.plot(q95,color='dimgrey',linestyle='solid') - ax.fill_between(np.arange(0,len(q5),1),list(q5),list(q95),facecolor='dimgrey',zorder=0, - linewidth=0,label='parameter uncertainty') - ax.plot(spot_setup.evaluation(),'r.',label='data') - ax.set_ylim(-50,450) - ax.set_xlim(0,729) - ax.legend() - fig.savefig('python_hymod.png',dpi=300) - -![Posterior model uncertainty](../img/python_hymod_simulation.png) -*Figure 1: Posterior model uncertainty of HYMOD.* - -## Plot convergence diagnostic -If you want to check the convergence of the DREAM algorithm: - - spotpy.analyser.plot_gelman_rubin(results, r_hat, fig_name = 'python_hymod_convergence.png') - -![Convergence diagnostic](../img/python_hymod_convergence.png) - -*Figure 2: Gelman-Rubin onvergence diagnostic of DREAM results.* - - -## Plot parameter uncertainty -Or if you want to check the posterior parameter distribution: - - parameters = spotpy.parameter.get_parameters_array(spot_setup) - - fig, ax = plt.subplots(nrows=5, ncols=2) - for par_id in range(len(parameters)): - plot_parameter_trace(ax[par_id][0], results, parameters[par_id]) - plot_posterior_parameter_histogram(ax[par_id][1], results, parameters[par_id]) - - ax[-1][0].set_xlabel('Iterations') - ax[-1][1].set_xlabel('Parameter range') - - plt.show() - fig.savefig('hymod_parameters.png',dpi=300) - -![Parameter distribution](../img/python_hymod_parameters.png) - -*Figure 3: Posterior parameter distribution of HYMOD.* diff --git a/docs/Calibration_with_SCE-UA.md b/docs/Calibration_with_SCE-UA.md new file mode 100644 index 00000000..fac6380e --- /dev/null +++ b/docs/Calibration_with_SCE-UA.md @@ -0,0 +1,88 @@ +# Calibration of HYMOD with SCE-UA + +This chapter shows you, how to calibrate an external hydrological model (HYMOD) with SPOTPY. + +We use the previously created hydrological model HYMOD spotpy_setup class as an [example](https://github.com/thouska/spotpy/blob/master/spotpy/examples/spot_setup_hymod_python.py), +to perform a calibration with the Shuffled Complex Evolution Algorithm (SCE-UA) algorithm. +For detailed information about the underlying theory, have a look at +Duan, Q., Sorooshian, S. and Gupta, V. K. (1994) +Optimal use of the SCE-UA global optimization method for calibrating watershed models, J. Hydrol.. + +First some relevant functions are imported Hymod example: + + import numpy as np + import spotpy + from spotpy.examples.spot_setup_hymod_python import spot_setup + import matplotlib.pyplot as plt + +As SCE-UA is minizing the objective function we need a objective function that gets better with decreasing values. +The Root Mean Squared Error objective function is a suitable choice: + + spot_setup=spot_setup(spotpy.objectivefunctions.rmse) + +Now we create a sampler, by using one of the implemented algorithms in SPOTPY. +The algorithm needs the inititalized spotpy setup class, and wants you to define a database name and database format. Here we create a SCEUA_hymod.csv file, +which will contain the returned likelihood function, the coresponding parameter set, the model results and a chain ID (SCE-UA is starting different independent +optimization complexes, which need to be in line, to receive robust results, they are herein called chains): + + sampler=spotpy.algorithms.sceua(spot_setup, dbname='SCEUA_hymod', dbformat='csv') + +To actually start the algorithm, spotpy needs some further details, like the maximum allowed number of repetitions, the number of chains used by dream (default = 5) and set the Gelman-Rubin convergence limit (default 1.2). +We further allow 100 runs after convergence is achieved: + + #Select number of maximum repetitions + rep=5000 + +We start the sampler and set some optional algorithm specific settings (check out the publication of SCE-UA for details): + + sampler.sample(rep, ngs=7, kstop=3, peps=0.1, pcento=0.1) + + +## Access the results +All gained results can be accessed from the SPOTPY csv-database: + + results = spotpy.analyser.load_csv_results('SCEUA_hymod') + +These results are stored in a simple structered numpy array, giving you all the flexibility at hand to analyse the results. +Herein we show a example how the objective function was minimized during sampling: + + fig= plt.figure(1,figsize=(9,5)) + plt.plot(results['like1']) + plt.show() + plt.ylabel('RMSE') + plt.xlabel('Iteration') + fig.savefig('SCEUA_objectivefunctiontrace.png',dpi=300) + +![Posterior model uncertainty](../img/SCEUA_objectivefunctiontrace.png) +*Figure 1: Root Mean Squared Error of Hymod model results during the optimization with SCE-UA algorithm.* + + +## Plot parameter uncertainty +Or if you want to check the best model run, we first need to find the run_id with the minimal objective function value + bestindex,bestobjf = spotpy.analyser.get_minlikeindex(results) + +Than we select the best model run: + + best_model_run = results[bestindex] + +And filter results for simulation results only: + + fields=[word for word in best_model_run.dtype.names if word.startswith('sim')] + best_simulation = list(best_model_run[fields]) + +The best model run can be plotted by using basic Matplotlib skills: + + fig= plt.figure(figsize=(16,9)) + ax = plt.subplot(1,1,1) + ax.plot(best_simulation,color='black',linestyle='solid', label='Best objf.='+str(bestobjf)) + ax.plot(spot_setup.evaluation(),'r.',markersize=3, label='Observation data') + plt.xlabel('Number of Observation Points') + plt.ylabel ('Discharge [l s-1]') + plt.legend(loc='upper right') + fig.savefig('SCEUA_best_modelrun.png',dpi=300) + +![SCEUA_best_modelrun](../img/SCEUA_best_modelrun.png) + +*Figure 2: Best model run of HYMOD, calibrated with SCE-UA, using RMSE as objective function.* + +The corresponding code is available for download [here](https://github.com/thouska/spotpy/blob/master/spotpy/examples/tutorial_sceua_hymod.py). \ No newline at end of file diff --git a/docs/How_to_link_a_model_to_SPOTPY.md b/docs/How_to_link_a_model_to_SPOTPY.md new file mode 100644 index 00000000..ccb08d07 --- /dev/null +++ b/docs/How_to_link_a_model_to_SPOTPY.md @@ -0,0 +1,126 @@ +## Linking a model with SPOTPY + +The SPOTPY package comes with several examples, which are designed to help you to set up your own research project. + +Herein we show, how to link an external model with SPOTPY. + +The hydrological model HYMOD serves as an example, which is a commonly used model in the hydrological community. +It takes meteorological data as input and produces discharge as output. + +Basically any model can be linked to SPOTPY as long it is somehow start-able with the Python programming language. +Linking a model to SPOTPY is done by following five consecutive steps, which are grouped in a spotpy_setup class. +If these steps are followed, SPOTPY can work with this class and analyse the model in an automatized way, giving you +powerful tools at hand. + +## Step 0: Import relevant packages and generate a spotpy setup class + + from spotpy.parameter import Uniform + from spotpy.objectivefunctions import rmse + from spotpy.examples.hymod_python.hymod import hymod + import os + +Generating a spotpy setup class is as easy as creating a class in Python: + + class spot_setup(object): + +## Step 1: Define the parameter of the model as class parameters + +Now, we can fill the setup class. In our case the model comes along with five parameters. +Here we use Python class parameter, to initialize the parameter for our model and make them readable for SPOTPY. +Needed information is the prior distribution (herein we choose a Uniform distribution), the minimum allowed parameter +setting (low) and the maximum allow setting (high). + + cmax = Uniform(low=1.0 , high=500) + bexp = Uniform(low=0.1 , high=2.0) + alpha = Uniform(low=0.1 , high=0.99) + Ks = Uniform(low=0.001 , high=0.10) + Kq = Uniform(low=0.1 , high=0.99) + +## Step 2: Write the def __init__ function, which takes care of any things which need to be done only once + +In this step, we can load all the data needed to start the model and set information, which might be needed in the following, +but want change during the model assessment. In our case, we want to set the path to the model: + + def __init__(self,obj_func=None): + #Find Path to Hymod on users system + self.owd = os.path.dirname(os.path.realpath(__file__)) + self.hymod_path = self.owd+os.sep+'hymod_python' + +Further we want to read in the forcing meteorological data and the observation data from file: + + self.PET,self.Precip = [], [] + self.date,self.trueObs = [], [] + #Load Observation data from file + climatefile = open(self.hymod_path+os.sep+'hymod_input.csv', 'r') + headerline = climatefile.readline()[:-1] + + #Read model forcing in working storage (this is done only ones) + if ';' in headerline: + self.delimiter = ';' + else: + self.delimiter = ',' + self.header = headerline.split(self.delimiter) + for line in climatefile: + values = line.strip().split(self.delimiter) + self.date.append(str(values[0])) + self.Precip.append(float(values[1])) + self.PET.append(float(values[2])) + self.trueObs.append(float(values[3])) + climatefile.close() + +Finaly, in this case a transformation factor is needed, which brings the output of the model from [mm/day] into [l s-1], where 1.783 is the catchment area [km²] + + self.Factor = 1.783 * 1000 * 1000 / (60 * 60 * 24) + +The keyword `obj_func` is not necessary in most cases. Herein we use it to change the function in Step 5, which makes this example flexible and applicable to various algorithms. + + self.obj_func = obj_func + +## Step 3: Write the def simulation function, which starts your model + +This function need to do three things. First it needs to be able to receive a parameter set from SPOTPY. +This will be an list of, in our case five values. Each value represents one setting for one parameter. +The HYMOD model receives these settings through a function. However, one could also write the settings into a file, which is read by the model (most common way). +Finally we start the model with the parameter set and collect the data, which are at the end returned. + + def simulation(self,x): + #Here the model is actualy started with a unique parameter combination that it gets from spotpy for each time the model is called + data = hymod(self.Precip, self.PET, x[0], x[1], x[2], x[3], x[4]) + sim=[] + for val in data: + sim.append(val*self.Factor) + #The first year of simulation data is ignored (warm-up) + return sim[366:] + +## Step 4: Write the def evaluation function, which returns your observation data + +Model assessment needs data, which tell you something about the reality that the model aims to reproduce. The HYMOD model produces discharge data, so herein we need to return observed discharge data + + def evaluation(self): + #The first year of simulation data is ignored (warm-up) + return self.trueObs[366:] + +## Step 5: Write the def objectivefunction, which returns how good the model fits the observation data + +In this last step, we compare the observed and the simulated data. So we need a function which can receive both information (simulation and evaluation). +SPOTPY will internally take care of times when this function needs to be called and will bring the corresponding data to this function. +Important is that the first keyword is handled as simulation/model results and the second as observation/evaluation results. +The SPOTPY package gives you several functions at hand, which compare the model and observation data with each other. +Please note that the selection of such a function highly influences the results gained with a optimization algorithm. Further some optimization algorithms +have specific needs at this function, as they wont to minimize or maximize the returned value by playing around with the parameter settings. Further some algorithms +need specific likelihood function. Herein we choose the Root Mean Squared Error (rmse) as an example. +A value of 0 would be a perfect fit of the simulation and observation data, +inf would be the worst possible fit. + + def objectivefunction(self,simulation,evaluation, params=None): + #SPOTPY expects to get one or multiple values back, + #that define the performance of the model run + if not self.obj_func: + # This is used if not overwritten by user + like = rmse(evaluation,simulation) + else: + #Way to ensure flexible spot setup class + like = self.obj_func(evaluation,simulation) + return like +Finally, save the file and check out the next example how to use the spotpy setup class with an algorithm. + +The shown example code is available [here](https://github.com/thouska/spotpy/blob/master/spotpy/examples/spot_setup_hymod_python.py) \ No newline at end of file diff --git a/docs/Sensitivity_analysis_with_FAST.md b/docs/Sensitivity_analysis_with_FAST.md new file mode 100644 index 00000000..72d2b874 --- /dev/null +++ b/docs/Sensitivity_analysis_with_FAST.md @@ -0,0 +1,63 @@ +# Sensitivity analysis of HYMOD with FAST + +SPOTPY gives you the opportunity to start a sensitivity analysis of your model. In this case, we included a global sensitivity analysis called "FAST" based on +Saltelli et al. (1999). + +The algorithm will tell you, how sensitive your parameters are based on whatever is given back by your objective function. Before you start to sample, you should know how how many +iterations you need to get an reliable information about your parameter. The number of iteration can be calculate after [Henkel et al. GLOBAL SENSITIVITY ANALYSIS OF NONLINEAR MATHEMATICAL MODELS - AN +IMPLEMENTATION OF TWO COMPLEMENTING VARIANCE-BASED ALGORITHMS, 2012] (https://www.informs-sim.org/wsc12papers/includes/files/con308.pdf): + +$$N = (1+4M^2(1+(k-2)d))k$$ + +with N = needed parameter iterations, M= inference factor (SPOTPY default M=4) and d= frequency step (SPOTPY default d=2) and k as the number of parameters of your model. + +In our case, we provide the hymod model as an [example](https://github.com/thouska/spotpy/blob/master/spotpy/examples/spot_setup_hymod_python.py) +You can start the simulation with + + +So the first step is as always to initialize the spotpy setup class example: + import spotpy + from spotpy.examples.spot_setup_hymod_python import spot_setup + spot_setup=spot_setup() + +As a next step we apply the above formular to determine the number of repetitions needed for the FAST analysis. In our case we have k=5 parameter +in the hymod model, so we need N=1345 repetitions to get reliable results: + + rep = 1345 + +And that's already it. Now we can start the sensitivity analysis: + sampler = spotpy.algorithms.fast(spot_setup, dbname='FAST_hymod', dbformat='csv') + sampler.sample(rep) + +This will take some time. Meanwhile SPOTPY will report you about the progress and the approximate duration. +Further SPOTPY will create a database for you with the name 'Fast_hymod.csv', depending on the users selection above. +This file will contain the return objectivefunction values, the sampled parameter sets, the simulation results and some information +about a chain (not relevant for FAST). + +These results can be loaded, e.g. directly from the sampler: + + results = sampler.get_data() + +OR, by loading the results from the database: + + results = spotpy.analyser.load_csv_results('FAST_hymod') + +Finally SPOTPY gives you some tools at hand, to analyse these results. E.g. if you want to determine, which are the three most important parameter of your +model, you can use the following: + + spotpy.analyser.plot_fast_sensitivity(results, number_of_sensitiv_pars=3) + + + analyser.plot_fast_sensitivity(results,number_of_sensitiv_pars=5) + +This will show you a Plot with the total Sensitivity index of all your parameters and in this case the five most sensitive parameters (can be adjusted). +Herein we use a already create spotpy setup class from the tutorial. The code for this class +is available + + +![FAST sensitivity results](../img/FAST_sensitivity.png) + +*Figure 1: Example output of a sensitivity analysis using the FAST (Fourier Amplitude Sensitivity Test). +Total sensitivity index is plotted for every model parameter. Sensitive parameters are plotted in blue, insensitive parameter in orange.* + +The corresponding code is available for download [here](https://github.com/thouska/spotpy/blob/master/spotpy/examples/tutorial_fast_hymod.py). \ No newline at end of file diff --git a/docs/index.md b/docs/index.md index 2afedac2..be990ea0 100644 --- a/docs/index.md +++ b/docs/index.md @@ -157,8 +157,10 @@ The results can be analysed with some pre-build statistics and plotting features rope.py # RObust Parameter Estimation fast.py # Fourier Amplitude Sensitivity Testing abc.py # Artificial Bee Colony - fscabc.py # Fitness Scaled Chaotic Artificial Bee Colony - dream.py # Differential Evolution Adaptive Metropolis + fscabc.py # Fitness Scaled Chaotic Artificial Bee Colony + dream.py # Differential Evolution Adaptive Metropolis + dds.py # Dynamically Dimensioned Search + parallel/ mpi.py #Basic Parralel Computing features diff --git a/mkdocs.yml b/mkdocs.yml index a68e12e2..5348b679 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -20,7 +20,10 @@ nav: - Rosenbrock: Rosenbrock.md - Griewank: Griewank.md - Ackley: Ackley.md - - Calibrating a hydrological model with DREAM: Calibrating a hydrological model with DREAM.md + - How to link a model to SPOTPY: How_to_link_a_model_to_SPOTPY.md + - Sensitivity analysis with FAST: Sensitivity_analysis_with_FAST.md + - Calibrating with SCE-UA: Calibration_with_SCE-UA.md + - Bayesian uncertainty analysis with DREAM: Bayesian_uncertainty_analysis_with_DREAM.md - Advanced_hints: Advanced_hints.md - Algorithm_guide: Algorithm_guide.md diff --git a/spotpy/algorithms/sceua.py b/spotpy/algorithms/sceua.py index 86679eb1..fe5e59bf 100644 --- a/spotpy/algorithms/sceua.py +++ b/spotpy/algorithms/sceua.py @@ -15,7 +15,7 @@ class sceua(_algorithm): """ - This class holds the Shuffled Complex Evolution Algortithm (SCE-UA) algorithm, + This class holds the Shuffled Complex Evolution Algorithm (SCE-UA) algorithm, based on: Duan, Q., Sorooshian, S. and Gupta, V. K. (1994) Optimal use of the SCE-UA global optimization method for calibrating watershed models, J. Hydrol. diff --git a/spotpy/analyser.py b/spotpy/analyser.py index 1c7db8de..f7eb0aa6 100644 --- a/spotpy/analyser.py +++ b/spotpy/analyser.py @@ -187,7 +187,7 @@ def get_minlikeindex(results): textv=text+str(index[0][0])+text2+value print(textv) - return index, minimum + return index[0][0], minimum def get_percentiles(results,sim_number=''): diff --git a/spotpy/examples/spot_setup_hymod_python.py b/spotpy/examples/spot_setup_hymod_python.py index 876f1f41..38f6d6a1 100644 --- a/spotpy/examples/spot_setup_hymod_python.py +++ b/spotpy/examples/spot_setup_hymod_python.py @@ -68,7 +68,7 @@ def evaluation(self): def objectivefunction(self,simulation,evaluation, params=None): #SPOTPY expects to get one or multiple values back, - #that define the performence of the model run + #that define the performance of the model run if not self.obj_func: # This is used if not overwritten by user like = rmse(evaluation,simulation) diff --git a/spotpy/examples/tutorial_fast_hymod.py b/spotpy/examples/tutorial_fast_hymod.py index 4ea0f05e..8c350cf4 100644 --- a/spotpy/examples/tutorial_fast_hymod.py +++ b/spotpy/examples/tutorial_fast_hymod.py @@ -5,15 +5,10 @@ :author: Tobias Houska -This class holds example code how to use the dream algorithm +This class holds example code how to use the FAST algorithm ''' -try: - import spotpy -except ImportError: - import sys - sys.path.append(".") - import spotpy +import spotpy from spotpy.examples.spot_setup_hymod_python import spot_setup @@ -23,7 +18,7 @@ spot_setup = spot_setup() #Select number of maximum repetitions - rep = 1000 + rep = 1345 #Start a sensitivity analysis sampler = spotpy.algorithms.fast(spot_setup, dbname='FAST_hymod', dbformat='csv') diff --git a/spotpy/examples/tutorial_sceua_hymod.py b/spotpy/examples/tutorial_sceua_hymod.py index 3fd3626f..171ce0d9 100644 --- a/spotpy/examples/tutorial_sceua_hymod.py +++ b/spotpy/examples/tutorial_sceua_hymod.py @@ -24,11 +24,7 @@ spot_setup=spot_setup(spotpy.objectivefunctions.rmse) #Select number of maximum allowed repetitions - rep=1000 - filename = 'SCEUA_hymod' - # Create the SCE-UA sampler of spotpy, alt_objfun is set to None to force SPOTPY - # to jump into the def objectivefunction in the spot_setup class (default is - # spotpy.objectivefunctions.rmse) + rep=5000 sampler=spotpy.algorithms.sceua(spot_setup, dbname='SCEUA_hymod', dbformat='csv') #Start the sampler, one can specify ngs, kstop, peps and pcento id desired @@ -37,114 +33,32 @@ # Load the results gained with the sceua sampler, stored in SCEUA_hymod.csv results = spotpy.analyser.load_csv_results('SCEUA_hymod') - - print(len(results), 'runs were saved.') + #Plot how the objective function was minimized during sampling fig= plt.figure(1,figsize=(9,5)) plt.plot(results['like1']) plt.show() plt.ylabel('RMSE') plt.xlabel('Iteration') - fig.savefig('hymod_objectivefunction.png',dpi=300) - - # Example plot to show the parameter distribution ###### - fig= plt.figure(2,figsize=(9,9)) - normed_value = 1 - - plt.subplot(5,2,1) - x = results['parcmax'] - for i in range(int(max(results['chain'])-1)): - index=np.where(results['chain']==i+1) #Ignores burn-in chain - plt.plot(x[index],'.') - plt.ylabel('cmax') - plt.ylim(spot_setup.cmax.minbound, spot_setup.cmax.maxbound) - - - plt.subplot(5,2,2) - x = x[int(len(results)*0.9):] #choose the last 10% of the sample - hist, bins = np.histogram(x, bins=20, density=True) - widths = np.diff(bins) - hist *= normed_value - plt.bar(bins[:-1], hist, widths) - plt.ylabel('cmax') - plt.xlim(spot_setup.cmax.minbound, spot_setup.cmax.maxbound) - - - plt.subplot(5,2,3) - x = results['parbexp'] - for i in range(int(max(results['chain'])-1)): - index=np.where(results['chain']==i+1) - plt.plot(x[index],'.') - plt.ylabel('bexp') - plt.ylim(spot_setup.bexp.minbound, spot_setup.bexp.maxbound) - - plt.subplot(5,2,4) - x = x[int(len(results)*0.9):] - hist, bins = np.histogram(x, bins=20, density=True) - widths = np.diff(bins) - hist *= normed_value - plt.bar(bins[:-1], hist, widths) - plt.ylabel('bexp') - plt.xlim(spot_setup.bexp.minbound, spot_setup.bexp.maxbound) - - - - plt.subplot(5,2,5) - x = results['paralpha'] - for i in range(int(max(results['chain'])-1)): - index=np.where(results['chain']==i+1) - plt.plot(x[index],'.') - plt.ylabel('alpha') - plt.ylim(spot_setup.alpha.minbound, spot_setup.alpha.maxbound) + fig.savefig('SCEUA_objectivefunctiontrace.png',dpi=300) - - plt.subplot(5,2,6) - x = x[int(len(results)*0.9):] - hist, bins = np.histogram(x, bins=20, density=True) - widths = np.diff(bins) - hist *= normed_value - plt.bar(bins[:-1], hist, widths) - plt.ylabel('alpha') - plt.xlim(spot_setup.alpha.minbound, spot_setup.alpha.maxbound) - - - plt.subplot(5,2,7) - x = results['parKs'] - for i in range(int(max(results['chain'])-1)): - index=np.where(results['chain']==i+1) - plt.plot(x[index],'.') - plt.ylabel('Ks') - plt.ylim(spot_setup.Ks.minbound, spot_setup.Ks.maxbound) - - - plt.subplot(5,2,8) - x = x[int(len(results)*0.9):] + # Plot the best model run + #Find the run_id with the minimal objective function value + bestindex,bestobjf = spotpy.analyser.get_minlikeindex(results) - hist, bins = np.histogram(x, bins=20, density=True) - widths = np.diff(bins) - hist *= normed_value - plt.bar(bins[:-1], hist, widths) - plt.ylabel('Ks') - plt.xlim(spot_setup.Ks.minbound, spot_setup.Ks.maxbound) + # Select best model run + best_model_run = results[bestindex] - - plt.subplot(5,2,9) - x = results['parKq'] - for i in range(int(max(results['chain'])-1)): - index=np.where(results['chain']==i+1) - plt.plot(x[index],'.') - plt.ylabel('Kq') - plt.ylim(spot_setup.Kq.minbound, spot_setup.Kq.maxbound) - plt.xlabel('Iterations') - - plt.subplot(5,2,10) - x = x[int(len(results)*0.9):] - hist, bins = np.histogram(x, bins=20, density=True) - widths = np.diff(bins) - hist *= normed_value - plt.bar(bins[:-1], hist, widths) - plt.ylabel('Kq') - plt.xlabel('Parameter range') - plt.xlim(spot_setup.Kq.minbound, spot_setup.Kq.maxbound) - plt.show() - fig.savefig('hymod_parameters.png',dpi=300) \ No newline at end of file + #Filter results for simulation results + fields=[word for word in best_model_run.dtype.names if word.startswith('sim')] + best_simulation = list(best_model_run[fields]) + + fig= plt.figure(figsize=(16,9)) + ax = plt.subplot(1,1,1) + ax.plot(best_simulation,color='black',linestyle='solid', label='Best objf.='+str(bestobjf)) + ax.plot(spot_setup.evaluation(),'r.',markersize=3, label='Observation data') + plt.xlabel('Number of Observation Points') + plt.ylabel ('Discharge [l s-1]') + plt.legend(loc='upper right') + fig.savefig('SCEUA_best_modelrun.png',dpi=300) +