diff --git a/docs/Advanced_hints.md b/docs/Advanced_hints.md index a86e1166..e8f90c0d 100644 --- a/docs/Advanced_hints.md +++ b/docs/Advanced_hints.md @@ -1,4 +1,5 @@ + # Advanced settings This chapter will show you, how to get the full power out of SPOTPY: diff --git a/docs/Bayesian_uncertainty_analysis_with_DREAM.md b/docs/Bayesian_uncertainty_analysis_with_DREAM.md index c0a4d104..cbeabb2b 100644 --- a/docs/Bayesian_uncertainty_analysis_with_DREAM.md +++ b/docs/Bayesian_uncertainty_analysis_with_DREAM.md @@ -92,7 +92,7 @@ If you want to see the remaining posterior model uncertainty: ax.legend() fig.savefig('python_hymod.png',dpi=300) -![Posterior model uncertainty](../img/python_hymod_simulation.png) +![Posterior simulation uncertainty](../img/DREAM_simulation_uncertainty_Hymod.png) *Figure 1: Posterior model uncertainty of HYMOD.* ## Plot convergence diagnostic @@ -100,7 +100,7 @@ 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) +![Convergence diagnostic](../img/DREAM_r_hat.png) *Figure 2: Gelman-Rubin onvergence diagnostic of DREAM results.* @@ -121,8 +121,8 @@ Or if you want to check the posterior parameter distribution: plt.show() fig.savefig('hymod_parameters.png',dpi=300) -![Parameter distribution](../img/python_hymod_parameters.png) +![DREAM_parameter_uncertainty_Hymod](../img/DREAM_parameter_uncertainty_Hymod.png) -*Figure 3: Posterior parameter distribution of HYMOD.* +*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). \ No newline at end of file diff --git a/docs/Calibration_with_SCE-UA.md b/docs/Calibration_with_SCE-UA.md index fac6380e..d522f4e7 100644 --- a/docs/Calibration_with_SCE-UA.md +++ b/docs/Calibration_with_SCE-UA.md @@ -57,14 +57,15 @@ Herein we show a example how the objective function was minimized during samplin *Figure 1: Root Mean Squared Error of Hymod model results during the optimization with SCE-UA algorithm.* -## Plot parameter uncertainty +## 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')] diff --git a/docs/Sensitivity_analysis_with_FAST.md b/docs/Sensitivity_analysis_with_FAST.md index 941b2932..4015d351 100644 --- a/docs/Sensitivity_analysis_with_FAST.md +++ b/docs/Sensitivity_analysis_with_FAST.md @@ -1,4 +1,5 @@ + # 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 @@ -8,7 +9,7 @@ The algorithm will tell you, how sensitive your parameters are based on whatever 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$$ + $$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. diff --git a/docs/img/FAST_sensitivity.png b/docs/img/FAST_sensitivity.png index 0aa706b6..e0815f85 100644 Binary files a/docs/img/FAST_sensitivity.png and b/docs/img/FAST_sensitivity.png differ diff --git a/docs/img/SCEUA_best_modelrun.png b/docs/img/SCEUA_best_modelrun.png index bf928fde..f2c47c70 100644 Binary files a/docs/img/SCEUA_best_modelrun.png and b/docs/img/SCEUA_best_modelrun.png differ diff --git a/docs/img/SCEUA_objectivefunctiontrace.png b/docs/img/SCEUA_objectivefunctiontrace.png index 549cb422..b2fc9ed3 100644 Binary files a/docs/img/SCEUA_objectivefunctiontrace.png and b/docs/img/SCEUA_objectivefunctiontrace.png differ diff --git a/spotpy/analyser.py b/spotpy/analyser.py index f7eb0aa6..690d363e 100644 --- a/spotpy/analyser.py +++ b/spotpy/analyser.py @@ -283,13 +283,13 @@ def plot_parameter_trace(ax, results, parameter): #THis function plots the parameter setting for each run for i in range(int(max(results['chain']))): index=np.where(results['chain']==i) - ax.plot(results['par'+parameter['name']][index],'.') + ax.plot(results['par'+parameter['name']][index],'.', markersize=2) ax.set_ylabel(parameter['name']) ax.set_ylim(parameter['minbound'], parameter['maxbound']) def plot_posterior_parameter_histogram(ax, results, parameter): - #This functing is plotting the last 10% of the results - ax.hist(results['par'+parameter['name']][int(len(results)*0.9):], + #This functing is the last 100 runs + ax.hist(results['par'+parameter['name']][-100:], bins =np.linspace(parameter['minbound'],parameter['maxbound'],20)) ax.set_ylabel('Density') ax.set_xlim(parameter['minbound'], parameter['maxbound']) @@ -468,7 +468,7 @@ def plot_fast_sensitivity(results,like_index=1,number_of_sensitiv_pars=10,fig_na import matplotlib.pyplot as plt parnames=get_parameternames(results) - fig=plt.figure(figsize=(16,6)) + fig=plt.figure(figsize=(9,6)) ax = plt.subplot(1,1,1) Si = get_sensitivity_of_fast(results, like_index=like_index) @@ -530,7 +530,7 @@ def plot_fast_sensitivity(results,like_index=1,number_of_sensitiv_pars=10,fig_na ax.plot(np.arange(-1,len(parnames)+1,1),[threshold]*(len(parnames)+2),'r--') ax.set_xlim(-0.5,len(parnames)-0.5) plt.tight_layout() - fig.savefig(fig_name,dpi=300) + fig.savefig(fig_name,dpi=150) def plot_heatmap_griewank(results,algorithms, fig_name='heatmap_griewank.png'): @@ -989,7 +989,7 @@ def plot_gelman_rubin(results, r_hat_values,fig_name='gelman_rub.png'): Output: Plot as seen for e.g. in (Sadegh and Vrugt 2014)''' import matplotlib.pyplot as plt - fig= plt.figure(figsize=(12,16)) + fig= plt.figure(figsize=(9,6)) ax1 = plt.subplot(2,1,1) for i in range(int(max(results['chain']))+1): index=np.where(results['chain']==i) @@ -1008,7 +1008,7 @@ def plot_gelman_rubin(results, r_hat_values,fig_name='gelman_rub.png'): ax2.set_ylabel('R$^d$ - convergence diagnostic') ax2.set_xlabel('Number of chainruns') ax2.legend() - fig.savefig(fig_name,dpi=300) + fig.savefig(fig_name,dpi=150) def gelman_rubin(x): '''NOT USED YET''' diff --git a/spotpy/examples/tutorial_dream_hymod.py b/spotpy/examples/tutorial_dream_hymod.py index 36de14fc..57da1852 100644 --- a/spotpy/examples/tutorial_dream_hymod.py +++ b/spotpy/examples/tutorial_dream_hymod.py @@ -40,8 +40,8 @@ acceptance_test_option = 6 sampler=spotpy.algorithms.dream(spot_setup, dbname='DREAM_hymod', dbformat='csv') - r_hat = sampler.sample(rep, nChains, nCr, eps, convergence_limit, - runs_after_convergence,acceptance_test_option) + #r_hat = sampler.sample(rep, nChains, nCr, eps, convergence_limit, + # runs_after_convergence,acceptance_test_option) @@ -53,7 +53,7 @@ # Example plot to show remaining parameter uncertainty # - fig= plt.figure(figsize=(16,9)) + fig= plt.figure(figsize=(9,6)) ax = plt.subplot(1,1,1) q5,q25,q75,q95=[],[],[],[] for field in fields: @@ -63,16 +63,16 @@ 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.plot(spot_setup.evaluation(),color='red', markersize=2,label='data') ax.set_ylim(-50,450) ax.set_xlim(0,729) ax.legend() - fig.savefig('python_hymod.png',dpi=300) + fig.savefig('DREAM_simulation_uncertainty_Hymod.png',dpi=150) ######################################################### # Example plot to show the convergence ################# - spotpy.analyser.plot_gelman_rubin(results, r_hat) + spotpy.analyser.plot_gelman_rubin(results, r_hat, fig_name='DREAM_r_hat.png') ######################################################## @@ -82,6 +82,8 @@ parameters = spotpy.parameter.get_parameters_array(spot_setup) fig, ax = plt.subplots(nrows=5, ncols=2) + fig.set_figheight(9) + fig.set_figwidth(9) 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]) @@ -90,5 +92,6 @@ ax[-1][1].set_xlabel('Parameter range') plt.show() - fig.savefig('hymod_parameters.png',dpi=300) + fig.tight_layout() + fig.savefig('DREAM_parameter_uncertainty_Hymod.png',dpi=300) ####################################################### \ No newline at end of file diff --git a/spotpy/examples/tutorial_sceua_hymod.py b/spotpy/examples/tutorial_sceua_hymod.py index 171ce0d9..abd2b380 100644 --- a/spotpy/examples/tutorial_sceua_hymod.py +++ b/spotpy/examples/tutorial_sceua_hymod.py @@ -35,12 +35,12 @@ #Plot how the objective function was minimized during sampling - fig= plt.figure(1,figsize=(9,5)) + fig= plt.figure(1,figsize=(9,6)) plt.plot(results['like1']) plt.show() plt.ylabel('RMSE') plt.xlabel('Iteration') - fig.savefig('SCEUA_objectivefunctiontrace.png',dpi=300) + fig.savefig('SCEUA_objectivefunctiontrace.png',dpi=150) # Plot the best model run #Find the run_id with the minimal objective function value @@ -53,12 +53,12 @@ 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)) + fig= plt.figure(figsize=(9,6)) 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) + fig.savefig('SCEUA_best_modelrun.png',dpi=150)