Fix figure plotting and printing in Tutorials

spotpy/algorithms/_algorithm.py

@@ -45,7 +45,7 @@ def __init__(self, repetitions, algorithm_name, optimization_direction, parnames
             print('The objective function will be minimized')
         if optimization_direction == 'maximize':
             self.compare = self.maximizer
-            print('The objective function will be minimized')
+            print('The objective function will be maximized')
         if optimization_direction == 'grid':
             self.compare = self.grid
 

spotpy/algorithms/padds.py

@@ -182,6 +182,7 @@ def sample(self, repetitions, trials=1, initial_objs=np.array([]), initial_param
         # Spotpy will not need this values.
         debug_results = []
         print('Starting the PADDS algotrithm with ' + str(repetitions) + ' repetitions...')
+        print('WARNING: THE PADDS algorithm as implemented in SPOTPY is in an beta stage and not ready for production use!')
         self.set_repetiton(repetitions)
         self.number_of_parameters = len(self.best_value.parameters) # number_of_parameters is the amount of parameters
 

spotpy/analyser.py

@@ -280,7 +280,7 @@ def get_posterior(results,percentage=10, maximize=True):
     return results[index]
 
 def plot_parameter_uncertainty(posterior_results,evaluation, fig_name='Posterior_parameter_uncertainty.png'):
-    import pylab as plt
+    import matplotlib.pyplot as plt
 
     simulation_fields = get_simulation_fields(posterior_results)
     fig= plt.figure(figsize=(16,9))
@@ -781,7 +781,7 @@ def plot_bestmodelrun(results,evaluation,fig_name ='Best_model_run.png'):
     Returns:
         figure. Plot of the simulation with the maximum objectivefunction value in the result array as a blue line and dots for the evaluation data.
     """
-    import pylab as plt
+    import matplotlib.pyplot as plt
     fig= plt.figure(figsize=(16,9))
     for i in range(len(evaluation)):
         if evaluation[i] == -9999:

spotpy/examples/MyOwnDatabase.txt

@@ -0,0 +1 @@
+obj_functions	parameters	simulations

spotpy/examples/hymod_exe/Param.in

@@ -1 +1 @@
-306.82853 1.20622 0.63259 0.05897 0.60575 
+401.08218 0.23071 0.43731 0.06884 0.16352 

spotpy/examples/spot_setup_hymod_python.py

@@ -12,31 +12,35 @@
 from __future__ import print_function
 from __future__ import unicode_literals
 
-import spotpy
+from spotpy.parameter import Uniform
+from spotpy.objectivefunctions import rmse
 from spotpy.examples.hymod_python.hymod import hymod
 import os
 
 class spot_setup(object):
-    cmax  = spotpy.parameter.Uniform(low=1.0 , high=500,  optguess=412.33)
-    bexp  = spotpy.parameter.Uniform(low=0.1 , high=2.0,  optguess=0.1725)
-    alpha = spotpy.parameter.Uniform(low=0.1 , high=0.99, optguess=0.8127)
-    Ks    = spotpy.parameter.Uniform(low=0.001 , high=0.10, optguess=0.0404)
-    Kq    = spotpy.parameter.Uniform(low=0.1 , high=0.99, optguess=0.5592)
+    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)
     #fake1 =spotpy.parameter.Uniform(low=0.1 , high=10, optguess=0.5592)
     #fake2 =spotpy.parameter.Uniform(low=0.1 , high=10, optguess=0.5592)
 
-    def __init__(self, _used_algorithm = 'default'):
-        self._used_algorithm = _used_algorithm  
+    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) 
         #Load Observation data from file
         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'
         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:
@@ -48,27 +52,27 @@ def __init__(self, _used_algorithm = 'default'):
             self.Precip.append(float(values[1]))
             self.PET.append(float(values[2]))
             self.trueObs.append(float(values[3]))
-
         climatefile.close()
 
-
     def simulation(self,x):
+        #Here the model is actualy startet with one paramter combination
         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:]
 
     def evaluation(self):
         return self.trueObs[366:]
 
     def objectivefunction(self,simulation,evaluation, params=None):
-        if self._used_algorithm == 'sceua':
-            like = spotpy.objectivefunctions.rmse(evaluation,simulation)
-        elif self._used_algorithm == 'dream':
-            like = spotpy.objectivefunctions.log_p(evaluation,simulation)
-        elif self._used_algorithm == 'default':
-            like = spotpy.objectivefunctions.nashsutcliffe(evaluation,simulation)
+        #SPOTPY expects to get one or multiple values back, 
+        #that define the performence of the model run
+        if not self.obj_func:
+            # This is used if not overwritten by user
+            like = rmse(evaluation,simulation)
         else:
-            like = spotpy.likelihoods.gaussianLikelihoodMeasErrorOut(evaluation,simulation)    
+            #Way to ensure on flexible spot setup class
+            like = self.obj_func(evaluation,simulation)    
         return like

spotpy/examples/spot_setup_rosenbrock.py

@@ -9,10 +9,10 @@
 from __future__ import division
 from __future__ import print_function
 from __future__ import unicode_literals
-import numpy as np
 
+import numpy as np
 from spotpy.parameter import Uniform
-from spotpy.objectivefunctions import rmse, log_p
+from spotpy.objectivefunctions import rmse
 
 class spot_setup(object):
     """
@@ -23,8 +23,10 @@ class spot_setup(object):
     x = Uniform(-10, 10, 1.5, 3.0, -10, 10, doc='x value of Rosenbrock function')
     y = Uniform(-10, 10, 1.5, 3.0, -10, 10, doc='y value of Rosenbrock function')
     z = Uniform(-10, 10, 1.5, 3.0, -10, 10, doc='z value of Rosenbrock function')
-    def __init__(self,used_algorithm='default'):
-        self.used_algorithm =used_algorithm
+
+    def __init__(self,obj_func=None):
+        self.obj_func = obj_func
+
     def simulation(self, vector):
         x=np.array(vector)
         simulations= [sum(100.0 * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0)]
@@ -35,10 +37,12 @@ def evaluation(self):
         return observations
 
     def objectivefunction(self, simulation, evaluation):
-        if self.used_algorithm == 'sceua' or self.used_algorithm == 'abc' or self.used_algorithm == 'fscabc':
-            objectivefunction = rmse(evaluation=evaluation, simulation=simulation)
-        elif self.used_algorithm == 'dream' or self.used_algorithm == 'demcz' or self.used_algorithm == 'mcmc':
-            objectivefunction = log_p(evaluation=evaluation, simulation=simulation)
+        #SPOTPY expects to get one or multiple values back, 
+        #that define the performence of the model run
+        if not self.obj_func:
+            # This is used if not overwritten by user
+            like = rmse(evaluation,simulation)
         else:
-            objectivefunction = - rmse(evaluation=evaluation, simulation=simulation)
-        return objectivefunction
+            #Way to ensure on flexible spot setup class
+            like = self.obj_func(evaluation,simulation)    
+        return like

spotpy/examples/tutorial_Parameterlist_iterator.py

@@ -15,11 +15,13 @@
 
 
 class spot_setup(object):
-    slow = 1000
+
     def __init__(self):
+        self.slow = 1000
         self.params = [spotpy.parameter.List('x',[1,2,3,4,5,6,7,8,9,0]), #Give possible x values as a List
                        spotpy.parameter.List('y',[0,1,2,5,6,7,8,9,0,1])  #Give possible y values as a List
-                       ]
+
+                   ]
     def parameters(self):
         return spotpy.parameter.generate(self.params)
 

spotpy/examples/tutorial_dds_hymod.py

@@ -17,7 +17,7 @@
     import spotpy
 
 from spotpy.examples.spot_setup_hymod_python import spot_setup
-import pylab as plt
+import matplotlib.pyplot as plt
 
 
 if __name__ == "__main__":
@@ -27,125 +27,73 @@
     # Initialize the Hymod example
     # In this case, we tell the setup which algorithm we want to use, so
     # we can use this exmaple for different algorithms
-    spot_setup=spot_setup()
+    spot_setup=spot_setup(users_objective_function=spotpy.objectivefunctions.nashsutcliffe)
 
     #Select number of maximum allowed repetitions
     rep=1000
 
     # 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) 
-    sampler=spotpy.algorithms.dds(spot_setup, dbname='DDS_hymod', dbformat='csv', alt_objfun=None)
+    sampler=spotpy.algorithms.dds(spot_setup, dbname='DDS_hymod', dbformat='csv')
 
     #Start the sampler, one can specify ngs, kstop, peps and pcento id desired
     sampler.sample(rep)
     results = sampler.getdata()
-    print(results)
+
     fig= plt.figure(1,figsize=(9,5))
     plt.plot(results['like1'])
     plt.show()
-    plt.ylabel('RMSE')
+    plt.ylabel('Objective function')
     plt.xlabel('Iteration')
     fig.savefig('hymod_objectivefunction.png',dpi=300)
 
     # Example plot to show the parameter distribution ###### 
+    def plot_parameter_trace(ax, results, parameter):
+        ax.plot(results['par'+parameter.name],'.')
+        ax.set_ylabel(parameter.name)
+        ax.set_ylim(parameter.minbound, parameter.maxbound)
+
+    def plot_parameter_histogram(ax, results, parameter):
+        #chooses the last 20% of the sample
+        ax.hist(results['par'+parameter.name][int(len(results)*0.9):], 
+                 bins =np.linspace(parameter.minbound,parameter.maxbound,20))
+        ax.set_ylabel('Density')
+        ax.set_xlim(parameter.minbound, parameter.maxbound)
+
     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)
+    ax1 = plt.subplot(5,2,1)
+    plot_parameter_trace(ax1, results, spot_setup.cmax)
 
+    ax2 = plt.subplot(5,2,2)
+    plot_parameter_histogram(ax2,results, spot_setup.cmax)
 
+    ax3 = plt.subplot(5,2,3)
+    plot_parameter_trace(ax3, results, spot_setup.bexp)
 
-    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)
-
-
-    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):]
+    ax4 = plt.subplot(5,2,4)
+    plot_parameter_histogram(ax4, results, spot_setup.bexp)
+
+    ax5 = plt.subplot(5,2,5)
+    plot_parameter_trace(ax5, results, spot_setup.alpha)
+
+    ax6 = plt.subplot(5,2,6)
+    plot_parameter_histogram(ax6, results, spot_setup.alpha)
 
-    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)
+    ax7 = plt.subplot(5,2,7)
+    plot_parameter_trace(ax7, results, spot_setup.Ks)
 
+    ax8 = plt.subplot(5,2,8)
+    plot_parameter_histogram(ax8, results, spot_setup.Ks)
+
+    ax9 = plt.subplot(5,2,9)
+    plot_parameter_trace(ax9, results, spot_setup.Kq)
+    ax9.set_xlabel('Iterations')
 
-    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')
+    ax10 = plt.subplot(5,2,10)
+    plot_parameter_histogram(ax10, results, spot_setup.Kq)
+    ax10.set_xlabel('Parameter range')
 
-    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)

spotpy/examples/tutorial_dream_hymod.py

@@ -12,14 +12,14 @@
 import spotpy
 #from spotpy.examples.spot_setup_hymod_exe import spot_setup
 from spotpy.examples.spot_setup_hymod_python import spot_setup
-import pylab as plt
+import matplotlib.pyplot as plt
 
 
 if __name__ == "__main__":
     parallel ='seq'
     # Initialize the Hymod example (will only work on Windows systems)
     #spot_setup=spot_setup(parallel=parallel)
-    spot_setup=spot_setup(_used_algorithm='dream')
+    spot_setup=spot_setup(obj_func = spotpy.likelihoods.gaussianLikelihoodMeasErrorOut)
 
     # Create the Dream sampler of spotpy, al_objfun is set to None to force SPOTPY
     # to jump into the def objectivefunction in the spot_setup class (default is

spotpy/examples/tutorial_fast_hymod.py

@@ -14,8 +14,7 @@
     import sys
     sys.path.append(".")
     import spotpy
-from spotpy.examples.spot_setup_hymod_exe import spot_setup
-#from spotpy.examples.spot_setup_hymod_python import spot_setup
+from spotpy.examples.spot_setup_hymod_python import spot_setup
 
 
 if __name__ == "__main__":

spotpy/examples/tutorial_likelihood.py

@@ -350,7 +350,7 @@ def objectivefunction(self, simulation=simulation, evaluation=evaluation, params
 sampler = spotpy.algorithms.lhs(spot_setup, dbname='RosenMC', dbformat='csv', alt_objfun=None)
 sampler.sample(rep)
 results.append(sampler.getdata())
-import pprint
+
 import pickle
 
 pickle.dump(results, open('mle_LimitsOfAcceptability_v2.p', 'wb'))