Package mystic :: Module scemtools

Module scemtools

source code

Implements the "Shuffled Complex Evolution Metropolis" Algoritm of Vrugt et al. [1]

Reference:

[1] Jasper A. Vrugt, Hoshin V. Gupta, Willem Bouten, and Soroosh Sorooshian A Shuffled Complex Evolution Metropolis algorithm for optimization and uncertainty assessment of hydrologic model parameters, WATER RESOURCES RESEARCH, VOL. 39, NO. 8, 1201, doi:10.1029/2002WR001642, 2003 http://www.agu.org/pubs/crossref/2003/2002WR001642.shtml

[2] Vrugt JA, Nuallain , Robinson BA, Bouten W, Dekker SC, Sloot PM Application of parallel computing to stochastic parameter estimation in environmental models, Computers & Geosciences, Vol. 32, No. 8. (October 2006), pp. 1139-1155. http://www.science.uva.nl/research/scs/papers/archive/Vrugt2006b.pdf

Functions
 
multinormal_pdf(mean, var)
var must be symmetric positive definite
source code
 
sequential_deal(inarray, n)
inarray: should be a set of N objects (the objects can be vectors themselves, but inarray should be index-able like a list.
source code
 
sort_and_deal(cards, target, nplayers) source code
 
sort_ab_with_b(a, b, ord=-1)
default is descending...
source code
 
sort_complex0(c, a) source code
 
sort_complex(c, a) source code
 
myinsert(a, x) source code
 
sort_complex2(c, a)
c and a are partially sorted (either the last one is bad, or the first one)
source code
 
update_complex(Ck, ak, c, a, pos)
ak is sorted (descending)
source code
 
scem(Ck, ak, Sk, Sak, target, cn)
This is the SCEM algorithm starting from line [35] of the reference [1].
source code
 
remix(Cs, As)
Mixing and dealing the complexes.
source code
Variables
  TYPE_DOUBLE = numpy.lib.getlimits.finfo(numpy.lib.getlimits.nu...
  TINY = TYPE_DOUBLE.tiny
Function Details

sequential_deal(inarray, n)

source code 
  • inarray: should be a set of N objects (the objects can be vectors themselves, but inarray should be index-able like a list. It is coerced into a numpy array because the last operations requires that it is also indexable by a 'list.'
  • it should have a length divisble by n, otherwise the reshape will fail (this is a feature !)
  • sequential_deal(range(20), 5) wil return a 5 element list, each element being a 4-list of index. (see below)
>>> for l in sequential_deal(range(20),5):
...    print l
...
[ 0  5 10 15]
[ 1  6 11 16]
[ 2  7 12 17]
[ 3  8 13 18]
[ 4  9 14 19]

scem(Ck, ak, Sk, Sak, target, cn)

source code 

This is the SCEM algorithm starting from line [35] of the reference [1].

- [inout] Ck is the kth 'complex' with m points. This should be an m by n array
  n being the dimensionality of the density function. i.e., the data
  are arranged in rows. 

  Ck is assumed to be sorted according to the target density.

- [inout] ak, the density of the points of Ck.

- [inout] Sk, the entire chain. (should be a list)

- [inout] Sak, the cost of the entire chain (should be a list)

  Sak would be more convenient to use if it is a numpy array, but we need
  to append to it frequently.

- [in] target: target density function

- [in] cn: jumprate. (see Paragraph 37 of [1.]

- The invariants: ak is always aligned with Ck, and are the cost of Ck
- Similarly, Sak is always aligned with Sk in the same way.

- On return... sort order in Ck/ak is destroyed. but see sort_complex2

    

remix(Cs, As)

source code 

Mixing and dealing the complexes. The types of Cs and As are very important....


Variables Details

TYPE_DOUBLE

Value:
numpy.lib.getlimits.finfo(numpy.lib.getlimits.numeric.double)