Does spotpy's Dream not exactly Vrugt, J. A. (2016) proposed Dream?

CHANGELOG.md

@@ -2,8 +2,10 @@
 
 All notable changes to this project will be documented in this file.
 
-## [Unreleased Changes](https://github.com/thouska/spotpy/compare/v1.5.13...master)
+## [Unreleased Changes](https://github.com/thouska/spotpy/compare/v1.5.15...master)
 
+## Spotpy Version [1.5.15] (https://github.com/thouska/spotpy/compare/v1.5.14..v.1.5.15) (2022-06-21)
+## Spotpy Version [1.5.14] (https://github.com/thouska/spotpy/compare/v1.5.13..v.1.5.14) (2020-10-09)
 ## Spotpy Version [1.5.13](https://github.com/thouska/spotpy/compare/v1.5.12...v1.5.13) (2020-09-07)
 
 * Introducing package dependencies as requested [#249](https://github.com/thouska/spotpy/issues/249)

setup.py

@@ -5,7 +5,7 @@
 
 setup(
   name = 'spotpy',
-  version = '1.5.14',
+  version = '1.5.15',
   description = 'A Statistical Parameter Optimization Tool',
   long_description=open(os.path.join(os.path.dirname(__file__),
                                        "README.rst")).read(),
@@ -27,5 +27,6 @@
         'Programming Language :: Python',
         'Programming Language :: Python :: 3.7',
         'Programming Language :: Python :: 3.8',
+        'Programming Language :: Python :: 3.9',
         'Topic :: Software Development :: Libraries :: Python Modules'],
         )

spotpy/__init__.py

@@ -41,4 +41,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.14'
+__version__ = '1.5.15'

spotpy/algorithms/dream.py

@@ -100,37 +100,59 @@ def check_par_validity_reflect(self, par):
             print('ERROR: Bounds have not the same lenghts as Parameterarray')
         return par
 
-    def _get_gamma(self,N):
+    def _get_gamma(self,newN, nchain_pairs):
         #N = Number of parameters
         p = np.random.uniform(low=0,high=1)
-        if p >=0.2:
-            gamma = 2.38/np.sqrt(2*int(N))#/self.gammalevel
+        if p >= 0.2:
+            # d_star is the dimension of subspace of parameters to jump
+            d_star = newN.count(True)
+            gamma = 2.38/np.sqrt(2*nchain_pairs*d_star)  # /self.gammalevel
         else:
             gamma = 1
         return gamma
 
-    def get_other_random_chains(self,cur_chain):
-        valid=False        
-        while valid == False:         
-            random_chain1 = np.random.randint(0,self.nChains)
-            random_chain2 = np.random.randint(0,self.nChains)
-            if random_chain1!=cur_chain and random_chain2!=cur_chain and random_chain1!=random_chain2:
-                valid=True
-        return random_chain1, random_chain2
+    def get_other_random_chains(self,cur_chain, nchain_pairs):
+        chain_pairs = []
+        selectable_chain = list(range(self.nChains))
+        selectable_chain.remove(cur_chain)
+
+        for i in range(nchain_pairs):
+            pair_ith = random.sample(selectable_chain, 2)
+            chain_pairs.append(pair_ith)
+            for chain in pair_ith:
+                selectable_chain.remove(chain)
+
+        return chain_pairs
 
-    def get_new_proposal_vector(self,cur_chain,newN,nrN):
-        gamma = self._get_gamma(nrN)
-        random_chain1,random_chain2 = self.get_other_random_chains(cur_chain)
+    def get_new_proposal_vector(self,cur_chain,newN,c):
+        nchain_pairs = random.randint(1, self.delta)
+        gamma = self._get_gamma(newN,nchain_pairs)
+        chain_pairs = self.get_other_random_chains(cur_chain, nchain_pairs)
         new_parameterset=[]        
         #position = self.chain_samples-1#self.nChains*self.chain_samples+self.chain_samples+cur_chain-1
         cur_par_set = list(self.bestpar[cur_chain][self.nChainruns[cur_chain]-1])
-        random_par_set1 = list(self.bestpar[random_chain1][self.nChainruns[random_chain1]-1])
-        random_par_set2 = list(self.bestpar[random_chain2][self.nChainruns[random_chain2]-1])
+        random_par_sets1 = []  # contain all random_par_set1
+        random_par_sets2 = []  # contain all random_par_set2
+
+        for i in range(nchain_pairs):
+            random_chain1 = chain_pairs[i][0]
+            random_chain2 = chain_pairs[i][1]
+            random_par_set1 = list(
+                self.bestpar[random_chain1][self.nChainruns[random_chain1]-1])
+            random_par_set2 = list(
+                self.bestpar[random_chain2][self.nChainruns[random_chain2]-1])
+            random_par_sets1.append(random_par_set1)
+            random_par_sets2.append(random_par_set2)
+
+        random_par_set1 = [sum(i) for i in zip(*random_par_sets1)]  # sum all random_par_set1
+        random_par_set2 = [sum(i) for i in zip(*random_par_sets2)]  # sum all random_par_set2
 
         for i in range(self.N):#Go through parameters
 
             if newN[i] == True:
-                new_parameterset.append(cur_par_set[i] + gamma*np.array(random_par_set1[i]-random_par_set2[i]) + np.random.normal(0,self.eps))
+                lambda_ =  np.random.uniform(-c,c)
+                new_parameterset.append(cur_par_set[i] + (1.0+lambda_)*gamma*np.array(
+                    random_par_set1[i]-random_par_set2[i]) + np.random.normal(0,self.eps))
             else:
                 new_parameterset.append(cur_par_set[i])
 
@@ -218,16 +240,17 @@ def get_r_hat(self, parameter_array):
 #            MR_stat = np.sqrt((n + 1) / n * R + (d - 1) / d)
             return R_stat#[R_stat, MR_stat]
 
-    def sample(self, repetitions,nChains=5, nCr=3, eps=10e-6, convergence_limit=1.2, runs_after_convergence=100,acceptance_test_option=6):
+    def sample(self, repetitions,nChains=5, nCr=3, delta=3, c=0.1, eps=10e-6, convergence_limit=1.2, runs_after_convergence=100,acceptance_test_option=6):
         self.set_repetiton(repetitions)
         print('Starting the DREAM algotrithm with '+str(repetitions)+ ' repetitions...')
-        if nChains <3:
-            print('Please use at least n=3 chains!')
+        if nChains <2*delta+1:
+            print('Please use at least n=2*delta+1 chains!')
             return None
         # Prepare storing MCMC chain as array of arrays.
         # define stepsize of MCMC.
         self.repetitions = int(repetitions)
         self.nChains = int(nChains)
+        self.delta = delta
         #Ensure initialisation of chains and database
         self.burnIn = self.nChains
         self.stepsizes = self.parameter()['step']  # array of stepsizes
@@ -272,7 +295,7 @@ def sample(self, repetitions,nChains=5, nCr=3, eps=10e-6, convergence_limit=1.2,
         nrN=1
         newN = [True]*self.N
         while self.iter < self.repetitions:
-            param_generator = ((curChain,self.get_new_proposal_vector(curChain,newN,nrN)) for curChain in range(int(self.nChains)))                
+            param_generator = ((curChain,self.get_new_proposal_vector(curChain,newN,c)) for curChain in range(int(self.nChains)))                
             for cChain,par,sim in self.repeat(param_generator):
                 pCr = np.random.randint(0,nCr)
                 ids=[]         

spotpy/examples/tutorial_dream_hymod.py

@@ -29,18 +29,20 @@
     #Select number of maximum repetitions
     rep=5000
 
-    # Select five chains and set the Gelman-Rubin convergence limit
-    nChains                = 4
+    # Select seven chains and set the Gelman-Rubin convergence limit
+    delta                  = 3
+    nChains                = 7
     convergence_limit      = 1.2
 
     # Other possible settings to modify the DREAM algorithm, for details see Vrugt (2016)
+    c                      = 0.1
     nCr                    = 3
     eps                    = 10e-6
     runs_after_convergence = 100
     acceptance_test_option = 6
 
     sampler=spotpy.algorithms.dream(spot_setup, dbname='DREAM_hymod', dbformat='csv', parallel=parallel)
-    r_hat = sampler.sample(rep, nChains, nCr, eps, convergence_limit, 
+    r_hat = sampler.sample(rep, nChains, nCr, delta, c, eps, convergence_limit, 
                            runs_after_convergence,acceptance_test_option)