Merge pull request #252 from lu-kas/master

Sample a predefined number of generations in SCEUA

setup.py

@@ -5,7 +5,7 @@
 
 setup(
   name = 'spotpy',
-  version = '1.5.10',
+  version = '1.5.11',
   description = 'A Statistical Parameter Optimization Tool',
   long_description=open(os.path.join(os.path.dirname(__file__),
                                        "README.rst")).read(),

spotpy/__init__.py

@@ -40,4 +40,4 @@
 from . import describe            # Contains some helper functions to describe samplers and set-ups
 from .hydrology import signatures # Quantifies goodness of fit between simulation and evaluation data with hydrological signatures
 
-__version__ = '1.5.10'
+__version__ = '1.5.11'

spotpy/algorithms/sceua.py

@@ -63,7 +63,6 @@ def __init__(self, *args, **kwargs):
         kwargs['optimization_direction'] = 'minimize'
         kwargs['algorithm_name'] = 'Shuffled Complex Evolution (SCE-UA) algorithm'
         super(sceua, self).__init__(*args, **kwargs)
-
 
     def simulate(self, id_params_tuple):
         """This overwrites the simple wrapper function of _algorithms.py
@@ -128,7 +127,7 @@ def simulate(self, id_params_tuple):
             # Replace the complex back into the population;
             return igs, likes, pars, sims, cx, cf, k1, k2, discarded_runs
 
-    def sample(self, repetitions, ngs=20, kstop=100, pcento=0.0000001, peps=0.0000001):
+    def sample(self, repetitions, ngs=20, kstop=100, pcento=0.0000001, peps=0.0000001, max_loop_inc=None):
         """
         Samples from parameter distributions using SCE-UA (Duan, 2004), 
         converted to python by Van Hoey (2011), restructured and parallelized by Houska et al (2015).
@@ -146,6 +145,8 @@ def sample(self, repetitions, ngs=20, kstop=100, pcento=0.0000001, peps=0.000000
             the percentage change allowed in the past kstop loops below which convergence is assumed to be achieved.
         peps: float
             Value of the normalized geometric range of the parameters in the population below which convergence is deemed achieved.
+        max_loop_inc: int
+            Number of loops executed at max in this function call
         """
         self.set_repetiton(repetitions)
         # Initialize SCE parameters:
@@ -163,6 +164,10 @@ def sample(self, repetitions, ngs=20, kstop=100, pcento=0.0000001, peps=0.000000
         bound = self.bu - self.bl  # np.array
         self.stochastic_parameters = bound != 0
         proceed = True
+
+        # burnin_only, needed to indictat if only the burnin phase should be run
+        burnin_only = False
+
         if self.breakpoint == 'read' or self.breakpoint == 'readandwrite':
             data_frombreak = self.read_breakdata(self.dbname)
             icall = data_frombreak[0]
@@ -177,7 +182,7 @@ def sample(self, repetitions, ngs=20, kstop=100, pcento=0.0000001, peps=0.000000
             icall = 0
             xf = np.zeros(npt)
 
-            print ('Starting burn-in sampling...')
+            print('Starting burn-in sampling...')
 
             # Burn in
             param_generator = ((rep, x[rep]) for rep in range(int(npt)))
@@ -194,6 +199,12 @@ def sample(self, repetitions, ngs=20, kstop=100, pcento=0.0000001, peps=0.000000
             idx = np.argsort(xf)
             xf = np.sort(xf)
             x = x[idx, :]
+
+            if max_loop_inc == 1:
+                burnin_only = True
+
+            print('Burn-in sampling completed...')
+
         else:
             raise ValueError("Don't know the breakpoint keyword {}".format(self.breakpoint))
 
@@ -221,16 +232,29 @@ def sample(self, repetitions, ngs=20, kstop=100, pcento=0.0000001, peps=0.000000
             print(
                 'THE POPULATION HAS CONVERGED TO A PRESPECIFIED SMALL PARAMETER SPACE')
 
-        print ('Burn-in sampling completed...')
+
+
 
         # Begin evolution loops:
         nloop = 0
         criter = []
         criter_change_pcent = 1e+5
-
-        self.repeat.setphase('ComplexEvo')
-        print ('Starting Complex Evolution...')
         proceed = True
+
+        # if only burnin, stop the following while loop to be started
+        # write brakpoint if only a single generation shall be computed and
+        # the main loop will not be executed
+        if burnin_only:
+            if self.breakpoint == 'write' or self.breakpoint == 'readandwrite':
+                work = (self.status.rep, (x, xf), gnrng)
+                self.write_breakdata(self.dbname, work)
+            proceed = False
+            print('ONLY THE BURNIN PHASE WAS COMPUTED')
+
+        else:
+            self.repeat.setphase('ComplexEvo')
+            print('Starting Complex Evolution...')
+
         while icall < repetitions and gnrng > peps and criter_change_pcent > pcento and proceed == True:
             nloop += 1
             print ('ComplexEvo loop #%d in progress...' % nloop)
@@ -260,10 +284,10 @@ def sample(self, repetitions, ngs=20, kstop=100, pcento=0.0000001, peps=0.000000
                         self.discarded_runs+=1
                         print('Skipping saving')
 
-                if self.breakpoint == 'write' or self.breakpoint == 'readandwrite'\
-                  and self.status.rep >= self.backup_every_rep:
-                    work = (self.status.rep, (x, xf), gnrng)
-                    self.write_breakdata(self.dbname, work)
+            if self.breakpoint == 'write' or self.breakpoint == 'readandwrite'\
+              and self.status.rep >= self.backup_every_rep:
+                work = (self.status.rep, (x, xf), gnrng)
+                self.write_breakdata(self.dbname, work)
 
             # End of Loop on Complex Evolution;
 
@@ -315,7 +339,13 @@ def sample(self, repetitions, ngs=20, kstop=100, pcento=0.0000001, peps=0.000000
                         'CONVERGENCY HAS ACHIEVED BASED ON OBJECTIVE FUNCTION CRITERIA!!!')
             elif self.status.stop:
                 proceed = False
-                break        
+                break
+
+            # stop, if max number of loop iteration was reached
+            elif max_loop_inc and nloop >= max_loop_inc:
+                proceed = False
+                print('THE MAXIMAL NUMBER OF LOOPS PER EXECUTION WAS REACHED')
+                break
 
         # End of the Outer Loops
         print('SEARCH WAS STOPPED AT TRIAL NUMBER: %d' % self.status.rep)
@@ -413,4 +443,4 @@ def _sampleinputmatrix(self, nrows, npars):
         x = np.zeros((nrows, npars))
         for i in range(nrows):
             x[i, :] = self.parameter()['random']
-        return x
+        return x

spotpy/examples/tutorial_sceua_hymod.py

@@ -21,18 +21,18 @@
     # 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(users_objective_function=spotpy.objectivefunctions.rmse)
+    spot_setup=spot_setup(spotpy.objectivefunctions.rmse)
 
     #Select number of maximum allowed repetitions
-    rep=500
+    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) 
     sampler=spotpy.algorithms.sceua(spot_setup, dbname='SCEUA_hymod', dbformat='csv')
 
     #Start the sampler, one can specify ngs, kstop, peps and pcento id desired
-    sampler.sample(rep)#,ngs=10, kstop=50, peps=0.1, pcento=0.1) 
+    sampler.sample(rep, ngs=7, kstop=3, peps=0.1, pcento=0.1) 
 
     # Load the results gained with the sceua sampler, stored in SCEUA_hymod.csv
     results = spotpy.analyser.load_csv_results('SCEUA_hymod')