Skip to content

Commit

Permalink
Merge branch 'master' into pr/246
Browse files Browse the repository at this point in the history
  • Loading branch information
thouska committed Oct 9, 2020
2 parents 4ae08b6 + afd996d commit f90a388
Show file tree
Hide file tree
Showing 39 changed files with 761 additions and 875 deletions.
28 changes: 24 additions & 4 deletions .travis.yml
Original file line number Diff line number Diff line change
@@ -1,10 +1,12 @@
language: python
python:
#- "2.7"
#- "3.5"
- "3.7"
# - "3.8" #Errors on Travis, build times out because no output was received...

# Uses cache to build faster
cache: pip

# this ubuntu distribution has sure the libopenmpi-dev packages available
dist: bionic

Expand All @@ -31,9 +33,27 @@ install:


script:
- py.test tests/test_* --cov spotpy --cov-report term-missing -v
- py.test tests/test_* --cov spotpy --cov-report term-missing -v
# - mpirun -c 4 python spotpy/examples/tutorial_parallel_computing_hymod.py 100
- pip uninstall spotpy -y

after_success:
- coveralls

jobs:
include:
- stage: Coveralls
name: Report Coveralls
after_success: python3 -m coveralls

# Publish spotpy on PyPi if taged
deploy:
provider: pypi
distributions: "sdist bdist_wheel"
on:
# In this case I want to deploy only when a tag is present...
tags: true
# ... and when tag is on master and respects the form "v0.0.0"
branch:
- master
- /v?(\d+\.)?(\d+\.)?(\*|\d+)$/
user: thouska
password : ${deploying}
10 changes: 10 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
# Changelog

All notable changes to this project will be documented in this file.

## [Unreleased Changes](https://github.com/thouska/spotpy/compare/v1.5.13...master)

## Spotpy Version [1.5.13](https://github.com/thouska/spotpy/compare/v1.5.12...v1.5.13) (2020-09-07)

* Introducing package dependencies as requested [#249](https://github.com/thouska/spotpy/issues/249)
* Introducing chanchelog.md as requested [#225](https://github.com/thouska/spotpy/issues/225)
25 changes: 2 additions & 23 deletions docs/Advanced_hints.md
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
<script type="text/javascript" src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML"></script>

# Advanced settings

This chapter will show you, how to get the full power out of SPOTPY:
Expand All @@ -20,7 +21,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.

Expand All @@ -38,29 +39,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:

Expand Down
128 changes: 128 additions & 0 deletions docs/Bayesian_uncertainty_analysis_with_DREAM.md
Original file line number Diff line number Diff line change
@@ -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 simulation uncertainty](../img/DREAM_simulation_uncertainty_Hymod.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/DREAM_r_hat.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)

![DREAM_parameter_uncertainty_Hymod](../img/DREAM_parameter_uncertainty_Hymod.png)

*Figure 3: Posterior parameter distribution of HYMOD. Plotting the last 100 repetitions of the algorithm as a histogram.*

The corresponding code is available for download [here](https://github.com/thouska/spotpy/blob/master/spotpy/examples/tutorial_dream_hymod.py).
89 changes: 89 additions & 0 deletions docs/Calibration_with_SCE-UA.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,89 @@
# 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 best model run
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).
Loading

0 comments on commit f90a388

Please sign in to comment.