Merge branch 'master' into New-Algorithm-DDS

* master: (22 commits)
  Added missing lines to allow for starting dream proposal vectors
  Update Version number upload to pypi
  Fix bug under mpi use
  Updates version number corresponds to upload on pypi
  Removed test for <Python3.6 due to deprectaed numpy version
  remove parameter interaction test for python 2
  Adopt test scripts to changes in examples
  Removed - sign from Rosenbrock example objectfivefunction
  Slight changes in sceua sampler and added corresponding tutorial
  Update _algorithm.py
  Work with None instead of np.NAN as this was not recognized
  Update __init__.py
  Update Version number, uploaded new pypi version
  Renamed keyword for saving switch
  Removes unfinished model runs from output file
  Enable automatic nan filtering for RMSE thouska#79
  Added comment
  Further version compability test
  Force pytest_cov down to v2.6
  Force decrease version of pytest_cov as v2.6.1 is deprecated
  ...

# Conflicts:
#	.travis.yml

.gitignore

@@ -9,6 +9,7 @@ __pycache__/
 #other files
 *.out
 *.in
+*.png
 site/*
 mkdocs.yml
 

.travis.yml

@@ -14,9 +14,8 @@ install:
   - pip install pandas==0.23.4
   - pip install scipy
   - pip install click
-  - pip install pytest==4.0.1
-  - pip install pytest-pep8
-  - pip install pytest-cov==2.6.0
+  # Need to force pytest-cov to v2.6 as current version (2.6.1) is deprecated
+  - pip install pytest pytest-pep8 pytest-cov==2.6
   # Use docutils to generate html describe
   - pip install docutils
   - pip install numba

setup.py

@@ -5,7 +5,7 @@
 
 setup(
   name = 'spotpy',
-  version = '1.3.30',
+  version = '1.4.2',
   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 .hydrology import signatures # Quantifies goodness of fit between simulation and evaluation data with hydrological signatures
 from . import tools
 
-__version__ = '1.3.30'
+__version__ = '1.4.2'

spotpy/algorithms/_algorithm.py

@@ -320,19 +320,20 @@ def update_params(self, params):
         return self.all_params
 
 
-    def postprocessing(self, rep, params, simulation, chains=1, save=True, negativlike=False): # TODO: rep not necessaray
+    def postprocessing(self, rep, params, simulation, chains=1, save_run=True, negativlike=False): # TODO: rep not necessaray
 
         params = self.update_params(params)
-        like = self.getfitness(simulation=simulation, params=params)
+        if negativlike is True:
+            like = -self.getfitness(simulation=simulation, params=params)
+        else:
+            like = self.getfitness(simulation=simulation, params=params)
+
         self.status(like, params)
         # Save everything in the database, if save is True
         # This is needed as some algorithms just want to know the fitness,
         # before they actually save the run in a database (e.g. sce-ua)
-        if save is True:
-            if negativlike is True:
-                self.save(-like, params, simulations=simulation, chains=chains)              
-            else:
-                self.save(like, params, simulations=simulation, chains=chains)    
+        if save_run is True and simulation is not None:
+            self.save(like, params, simulations=simulation, chains=chains)
         if type(like)==type([]):
             return like[0]
         else:        
@@ -378,8 +379,8 @@ def model_layer(q,all_params):
         sim_thread.join(self.sim_timeout)
 
         # If no result from the thread is given, i.e. the thread was killed from the watcher the default result is
-        # '[nan]' otherwise get the result from the thread
-        model_result = [np.NAN]
+        # '[nan]' and will not be saved. Otherwise get the result from the thread
+        model_result = None
         if not que.empty():
             model_result = que.get()
         return id, params, model_result

spotpy/algorithms/abc.py

@@ -103,7 +103,7 @@ def sample(self, repetitions, eb=48, a=(1 / 10), peps=0.0001, ownlimit=False, li
             work.append([like, randompar, like, randompar, c, p])
             icall +=1
             if self.status.stop:
-                print('Stopping samplig')
+                print('Stopping sampling')
                 break
 
         while icall < repetitions and gnrng > peps:
@@ -205,5 +205,6 @@ def sample(self, repetitions, eb=48, a=(1 / 10), peps=0.0001, ownlimit=False, li
 
             if gnrng < peps:
                 print(
-                    'THE POPULATION HAS CONVERGED TO A PRESPECIFIED SMALL PARAMETER SPACE')
+                    'THE POPULATION HAS CONVERGED TO A PRESPECIFIED SMALL PARAMETER SPACE AT RUN')
+                print(icall)
         self.final_call()

spotpy/algorithms/demcz.py

@@ -215,7 +215,9 @@ def sample(self, repetitions, nChains=3, burnIn=100, thin=1,
             # 3) and we have not done more than the maximum number of iterations
 
             while cur_iter < maxChainDraws:
+                print(cur_iter, burnIn)
                 if cur_iter == burnIn:
+                    print('starting')
                     history.start_sampling()
 
                 # every5th iteration allow a big jump
@@ -262,22 +264,15 @@ def sample(self, repetitions, nChains=3, burnIn=100, thin=1,
                     # each chain proposal
                     decisions, acceptance = self._metropolis_hastings(
                         currentLogPs, proposalLogPs, nChains)
-                    try:
-                        self._update_accepts_ratio(accepts_ratio_weighting, acceptance)
-                    except DEMCZError:
-                        pass
-
-
-
+                    self._update_accepts_ratio(accepts_ratio_weighting, acceptance)        
                     # choose from list of possible choices if 1d_decision is True at
                     # specific index, else use default choice
                     # np.choose(1d_decision[:,None], (list of possible choices, default
                     # choice)
                     save_likes=[]
                     save_pars=[]
                     save_sims=[]
-                    #print(len(self._logPs))
-
+
                     for curchain in range(nChains):
                         if decisions[curchain]:
                            save_likes.append(float(new_likelist[curchain]))
@@ -324,7 +319,8 @@ def sample(self, repetitions, nChains=3, burnIn=100, thin=1,
                         cur_iter = maxChainDraws
                         print(
                             'All chains fullfil the convergence criteria. Sampling stopped.')
-
+                cur_iter+=1
+
         # 3) finalize
         # only make the second half of draws available because that's the only
         # part used by the convergence diagnostic
@@ -417,6 +413,9 @@ def start_sampling(self):
     def sequence_histories(self):
         return self.group_sequence_histories('all')
 
+    def group_sequence_histories(self, name):
+        return self._sequence_histories[:, self.group_indicies[name], int(np.ceil(self.relevantHistoryStart)):self.relevantHistoryEnd]
+
     @property
     def nsequence_histories(self):
         return self.sequence_histories.shape[2]

spotpy/algorithms/sceua.py

@@ -65,13 +65,6 @@ def __init__(self, *args, **kwargs):
             kwargs['alt_objfun'] = 'rmse'
         super(sceua, self).__init__(*args, **kwargs)
 
-    def find_min_max(self):
-        randompar = self.parameter()['random']
-        for i in range(1000):
-            randompar = np.column_stack(
-                (randompar, self.parameter()['random']))
-        return np.amin(randompar, axis=1), np.amax(randompar, axis=1)
-
     def simulate(self, id_params_tuple):
         """This overwrites the simple wrapper function of _algorithms.py
         and makes a two phase mpi parallelization possbile:
@@ -84,7 +77,7 @@ def simulate(self, id_params_tuple):
 
         else:  # complex-evolution
             igs, x, xf, icall, cx, cf, sce_vars = id_params_tuple
-            self.npg, self.nopt, self.ngs, self.nspl, self.nps, self.bl, self.bu, self.status = sce_vars
+            self.npg, self.nopt, self.ngs, self.nspl, self.nps, self.bl, self.bu, self.status, self.stochastic_parameters = sce_vars
             # Partition the population into complexes (sub-populations);
             k1 = np.arange(self.npg, dtype=int)
             k2 = k1 * self.ngs + igs
@@ -192,11 +185,10 @@ def sample(self, repetitions, ngs=20, kstop=100, pcento=0.0000001, peps=0.000000
             param_generator = ((rep, x[rep]) for rep in range(int(npt)))
             for rep, randompar, simulations in self.repeat(param_generator):
                 # Calculate the objective function
-                like = self.postprocessing(icall, randompar, simulations, negativlike=True)              
+                like = self.postprocessing(icall, randompar, simulations,chains=0)              
                 xf[rep] = like
                 icall += 1
                 if self.status.stop:
-                    #icall = repetitions
                     print('Stopping samplig')
                     break
             # Sort the population in order of increasing function values;
@@ -250,15 +242,15 @@ def sample(self, repetitions, ngs=20, kstop=100, pcento=0.0000001, peps=0.000000
             cf = np.zeros((self.npg))
 
             sce_vars = [self.npg, self.nopt, self.ngs, self.nspl,
-                        self.nps, self.bl, self.bu, self.status]
+                        self.nps, self.bl, self.bu, self.status, self.stochastic_parameters]
             param_generator = ((rep, x, xf, icall, cx, cf, sce_vars)
                                for rep in range(int(self.ngs)))
             for igs, likes, pars, sims, cx, cf, k1, k2 in self.repeat(param_generator):
                 icall += len(likes)
                 x[k2, :] = cx[k1, :]
                 xf[k2] = cf[k1]
                 for i in range(len(likes)):
-                    like = self.postprocessing(icall+i, pars[i], sims[i], chains=i, negativlike=True)
+                    like = self.postprocessing(icall+i, pars[i], sims[i], chains=i+1)
                     if self.status.stop:
                         icall = repetitions
                         print('Stopping samplig')
@@ -381,7 +373,7 @@ def _cceua(self, s, sf, icall):
 
         ##    fnew = functn(self.nopt,snew);
         _, _, simulations = _algorithm.simulate(self, (icall, snew))
-        like = self.postprocessing(icall, snew, simulations, save=False)
+        like = self.postprocessing(icall, snew, simulations, save_run=False)
         fnew = like
         icall += 1
         if self.status.stop:
@@ -394,7 +386,7 @@ def _cceua(self, s, sf, icall):
             snew[constant_parameters] = sw[constant_parameters]
 
             _, _, simulations = _algorithm.simulate(self, (icall, snew))
-            like = self.postprocessing(icall, snew, simulations, save=False)
+            like = self.postprocessing(icall, snew, simulations, save_run=False)
             fnew = like
             icall += 1
             if self.status.stop:
@@ -406,7 +398,7 @@ def _cceua(self, s, sf, icall):
             if fnew > fw:
                 snew = self._sampleinputmatrix(1, self.nopt)[0]
                 _, _, simulations = _algorithm.simulate(self, (icall, snew))
-                like = self.postprocessing(icall, snew, simulations, save=False)
+                like = self.postprocessing(icall, snew, simulations, save_run=False)
                 fnew = like
                 icall += 1
 

spotpy/analyser.py

@@ -1016,7 +1016,7 @@ def plot_autocorellation(parameterdistribution,parametername):
     cnames=list(colors.cnames)
     fig=plt.figure(figsize=(16,9))
     ax = plt.subplot(1,1,1)
-    pd.tools.plotting.autocorrelation_plot(parameterdistribution)
+    pd.plotting.autocorrelation_plot(parameterdistribution)
     plt.savefig('Autocorellation'+str(parametername),dpi=300)
 
 

spotpy/database.py

@@ -112,7 +112,7 @@ def _make_header(self, simulations,parnames):
         #print(type(self.singular_data_lens[2]))        
         if self.save_sim:
             for i in range(len(simulations)):
-                if type(simulations[0]) == type([]) or type(simulations[0]) == type(np.array([])):
+                if isinstance(simulations[0], list) or type(simulations[0]) == type(np.array([])):
                     for j in range(len(simulations[i])):
                         self.header.extend(['simulation' + str(i+1)+'_'+str(j+1)])
                 else:

spotpy/examples/spot_setup_hymod_exe.py

@@ -98,5 +98,5 @@ def evaluation(self):
         return self.evals
 
     def objectivefunction(self,simulation,evaluation, params=None):
-        like = spotpy.objectivefunctions.nashsutcliffe(evaluation,simulation)     # Works good
+        like = spotpy.objectivefunctions.nashsutcliffe(evaluation,simulation)     # Works good for dream sampler
         return like

spotpy/examples/spot_setup_hymod_python.py

@@ -17,7 +17,16 @@
 import os
 
 class spot_setup(object):
-    def __init__(self):
+    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
@@ -41,16 +50,7 @@ def __init__(self):
             self.trueObs.append(float(values[3]))
 
         climatefile.close()
-        self.params = [spotpy.parameter.Uniform('cmax',low=1.0 , high=500,  optguess=412.33),
-                   spotpy.parameter.Uniform('bexp',low=0.1 , high=2.0,  optguess=0.1725),
-                   spotpy.parameter.Uniform('alpha',low=0.1 , high=0.99, optguess=0.8127),
-                   spotpy.parameter.Uniform('Ks',low=0.0 , high=0.10, optguess=0.0404),
-                   spotpy.parameter.Uniform('Kq',low=0.1 , high=0.99, optguess=0.5592),
-                   spotpy.parameter.Uniform('fake1',low=0.1 , high=10, optguess=0.5592),
-                   spotpy.parameter.Uniform('fake2',low=0.1 , high=10, optguess=0.5592)]
 
-    def parameters(self):
-        return spotpy.parameter.generate(self.params)
 
     def simulation(self,x):
         data = hymod(self.Precip, self.PET, x[0], x[1], x[2], x[3], x[4])
@@ -63,5 +63,8 @@ def evaluation(self):
         return self.trueObs[366:]
 
     def objectivefunction(self,simulation,evaluation, params=None):
-        like = spotpy.likelihoods.gaussianLikelihoodMeasErrorOut(evaluation,simulation)
+        if self._used_algorithm == 'sceua':
+            like = spotpy.objectivefunctions.rmse(evaluation,simulation)
+        else:
+            like = spotpy.likelihoods.gaussianLikelihoodMeasErrorOut(evaluation,simulation)    
         return like

spotpy/examples/spot_setup_rosenbrock.py

@@ -34,5 +34,5 @@ def evaluation(self):
         return observations
 
     def objectivefunction(self, simulation, evaluation):
-        objectivefunction = -rmse(evaluation=evaluation, simulation=simulation)
+        objectivefunction = rmse(evaluation=evaluation, simulation=simulation)
         return objectivefunction