Skip to content

Commit

Permalink
Merge pull request #285 from BeingHapppy/master
Browse files Browse the repository at this point in the history
Does spotpy's Dream not exactly Vrugt, J. A. (2016) proposed Dream?
  • Loading branch information
thouska authored Jun 21, 2022
2 parents 92c3aad + 2b844a7 commit e77e52c
Show file tree
Hide file tree
Showing 5 changed files with 55 additions and 27 deletions.
4 changes: 3 additions & 1 deletion CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
3 changes: 2 additions & 1 deletion setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(),
Expand All @@ -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'],
)
2 changes: 1 addition & 1 deletion spotpy/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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'
65 changes: 44 additions & 21 deletions spotpy/algorithms/dream.py
Original file line number Diff line number Diff line change
Expand Up @@ -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])

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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=[]
Expand Down
8 changes: 5 additions & 3 deletions spotpy/examples/tutorial_dream_hymod.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)


Expand Down

0 comments on commit e77e52c

Please sign in to comment.