PADDS tidy up tutorials, tests

spotpy/algorithms/_algorithm.py

@@ -244,8 +244,6 @@ def __init__(self, spot_setup, dbname=None, dbformat=None, dbinit=True,
         else:
             raise ImportError("Unsupported Python version")
 
-        print(objfunc_signature)
-
         objfunc_dynamic_arguments = {}
         for osp in (objfunc_signature):
             if osp == "self":

spotpy/algorithms/padds.py

@@ -142,7 +142,7 @@ def sample(self, repetitions, trials=1, x_initial=np.array([]), metric="ones"):
         # every iteration a map of all relevant values is stored, only for debug purpose.
         # Spotpy will not need this values.
         debug_results = []
-
+        print('Starting the PADDS algotrithm with ' + str(repetitions) + ' repetitions...')
         self.set_repetiton(repetitions)
         self.number_of_parameters = len(self.status.params) # number_of_parameters is the amount of parameters
 
@@ -209,7 +209,6 @@ def calc_initial_pareto_front(self, its):
             self.obj_func_current = self.getfitness(simulation=model_simulations, params=self.parameter_current)
             # First value will be used to initialize the values
             if i == 0:
-                print(self.pareto_front[0], np.append(self.obj_func_current, self.parameter_current))
                 self.pareto_front = np.vstack([self.pareto_front[0], np.append(self.obj_func_current, self.parameter_current)])
             else:
                 (self.pareto_front, dominance_flag) = nd_check(self.pareto_front, self.obj_func_current, self.parameter_current)

spotpy/examples/spot_setup_hymod_python.py

@@ -64,12 +64,4 @@ def evaluation(self):
         return self.trueObs[366:]
 
     def objectivefunction(self,simulation,evaluation, params=None):
-
-        #return -spotpy.likelihoods.gaussianLikelihoodMeasErrorOut(evaluation, simulation)
-        return [
-            #spotpy.likelihoods.gaussianLikelihoodMeasErrorOut(evaluation, simulation),
-            -spotpy.objectivefunctions.rmse(evaluation, simulation),
-            -spotpy.objectivefunctions.mse(evaluation, simulation),
-            #-spotpy.objectivefunctions.pbias(evaluation, simulation),
-            #spotpy.likelihoods.NashSutcliffeEfficiencyShapingFactor(evaluation, simulation)
-        ]
+        return -spotpy.likelihoods.gaussianLikelihoodMeasErrorOut(evaluation, simulation)

spotpy/examples/spot_setup_hymod_python_pareto.py

@@ -0,0 +1,73 @@
+'''
+Copyright 2015 by Tobias Houska
+This file is part of Statistical Parameter Estimation Tool (SPOTPY).
+
+:author: Tobias Houska
+
+This example implements the python version of hymod into SPOTPY.
+'''
+
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import spotpy
+from spotpy.examples.hymod_python.hymod import hymod
+import os
+import numpy as np
+
+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.0 , high=0.10, optguess=0.0404)
+    Kq    = spotpy.parameter.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  
+        #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 = [], []
+        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]
+
+        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()
+
+
+    def simulation(self,x):
+        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)
+        return sim[366:]
+
+    def evaluation(self):
+        return self.trueObs[366:]
+
+    def objectivefunction(self,simulation,evaluation, params=None):
+        return [
+            spotpy.likelihoods.gaussianLikelihoodMeasErrorOut(evaluation, simulation),
+            -spotpy.objectivefunctions.rmse(evaluation, simulation),
+            -spotpy.objectivefunctions.mse(evaluation, simulation),
+            -spotpy.objectivefunctions.pbias(evaluation, simulation),
+            #spotpy.likelihoods.NashSutcliffeEfficiencyShapingFactor(evaluation, simulation)
+        ]

spotpy/examples/tutorial_padds_hymod.py

@@ -16,7 +16,7 @@
     sys.path.append(".")
     import spotpy
 
-from spotpy.examples.spot_setup_hymod_python import spot_setup
+from spotpy.examples.spot_setup_hymod_python_pareto import spot_setup
 import pylab as plt
 
 
@@ -30,12 +30,12 @@
     spot_setup=spot_setup()
 
     #Select number of maximum allowed repetitions
-    rep=1000
+    rep=3000
 
     # 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.padds(spot_setup, dbname='padds_hymod', dbformat='csv', num_objs=3, alt_objfun=None)
+    sampler=spotpy.algorithms.padds(spot_setup, dbname='padds_hymod', dbformat='csv', alt_objfun=None)
 
     #Start the sampler, one can specify ngs, kstop, peps and pcento id desired
     print(sampler.sample(rep, metric="crowd_distance"))
@@ -45,8 +45,11 @@
 
 
     results = sampler.getdata()
+    from pprint import pprint
+    #pprint(results)
+    pprint(results['chain'])
 
-    for likno in range(1,4):
+    for likno in range(1,5):
         fig_like1 = plt.figure(1,figsize=(9,5))
         plt.plot(results['like'+str(likno)])
         plt.show()
@@ -101,6 +104,7 @@
 
     plt.subplot(5,2,5)
     x = results['paralpha']
+    print(x)
     for i in range(int(max(results['chain'])-1)):
         index=np.where(results['chain']==i+1)
         plt.plot(x[index],'.')