ENH: Add vector initializers

Add vector seeding for arrays with large states
Add implementation details to docstrings

Author: Sheppard, Kevin

randomstate/defaults.pxi

@@ -3,4 +3,5 @@ DEF RS_RNG_SEED = 1
 DEF RS_RNG_ADVANCEABLE = 0
 DEF RS_RNG_JUMPABLE = 0
 DEF RS_RNG_STATE_LEN = 4
-DEF RS_SEED_NBYTES = 2
+DEF RS_SEED_NBYTES = 2
+DEF RS_SEED_ARRAY_BITS = 64

randomstate/interface.pyx

@@ -13,14 +13,14 @@ cimport numpy as np
 cimport cython
 
 from libc cimport string
-from libc.stdint cimport (uint8_t, uint16_t, uint32_t, uint64_t, int8_t,
-                          int16_t, int32_t, int64_t, intptr_t)
+from libc.stdint cimport (uint8_t, uint16_t, uint32_t, uint64_t,
+                          int8_t, int16_t, int32_t, int64_t, intptr_t)
+
 from cpython cimport Py_INCREF
-from cpython.mem cimport PyMem_Malloc, PyMem_Free
 
 import randomstate
 from binomial cimport binomial_t
-from cython_overrides cimport PyFloat_AsDouble, PyInt_AsLong, PyErr_Occurred, PyErr_Clear
+from cython_overrides cimport PyFloat_AsDouble, PyInt_AsLong
 from randomstate.entropy import random_entropy
 
 np.import_array()
@@ -234,7 +234,9 @@ cdef class RandomState:
                     raise ValueError("Seed must be between 0 and 4294967295")
                 obj = obj.astype('L', casting='unsafe')
                 with self.lock:
-                    set_seed_by_array(&self.rng_state, <unsigned long *>np.PyArray_DATA(obj), np.PyArray_DIM(obj, 0))
+                    set_seed_by_array(&self.rng_state,
+                                      <unsigned long *>np.PyArray_DATA(obj),
+                                      np.PyArray_DIM(obj, 0))
             self._reset_state_variables()
 
     ELIF RS_RNG_SEED==1:
@@ -263,12 +265,36 @@ cdef class RandomState:
             --------
             RandomState
             """
-            if seed is not None:
-                if seed < 0:
-                    raise ValueError('seed < 0')
-            else:
-                self.__seed = seed = _generate_seed(RS_SEED_NBYTES)
-            set_seed(&self.rng_state, seed)
+            try:
+                if seed is not None:
+                    idx = operator.index(seed)
+                    if idx < 0:
+                        raise ValueError('seed < 0')
+                else:
+                    self.__seed = seed = _generate_seed(RS_SEED_NBYTES)
+                set_seed(&self.rng_state, seed)
+            except TypeError:
+                IF RS_SEED_ARRAY_BITS == 32:
+                    seed = np.asarray(seed).astype(np.int64, casting='safe')
+                    if ((seed > int(2**32 - 1)) | (seed < 0)).any():
+                        raise ValueError("Seed values must be between 0 and "
+                                         "4294967295  (2**32-1)")
+                    seed = seed.astype(np.uint32, casting='unsafe')
+                    with self.lock:
+                        set_seed_by_array(&self.rng_state,
+                                          <uint32_t *>np.PyArray_DATA(seed),
+                                          np.PyArray_DIM(seed, 0))
+                ELSE:
+                    seed = np.asarray(seed).astype(np.object, casting='safe')
+                    if ((seed > int(2**64 - 1)) | (seed < 0)).any():
+                        raise ValueError("Seed values must be between 0 and "
+                                         "18446744073709551616 (2**64-1)")
+                    seed = seed.astype(np.uint64, casting='unsafe')
+                    with self.lock:
+                        set_seed_by_array(&self.rng_state,
+                                          <uint64_t *>np.PyArray_DATA(seed),
+                                          np.PyArray_DIM(seed, 0))
+                self.__seed = seed
             self._reset_state_variables()
 
     ELSE:

randomstate/shims/dSFMT/dSFMT.pxi

@@ -4,6 +4,8 @@ DEF DSFMT_MEXP = 19937
 DEF DSFMT_N = 191 # ((DSFMT_MEXP - 128) / 104 + 1)
 DEF DSFMT_N_PLUS_1 = 192 # DSFMT_N + 1
 DEF RS_SEED_NBYTES = 1
+DEF RS_SEED_ARRAY_BITS = 32
+
 
 ctypedef uint32_t rng_state_t
 
@@ -103,11 +105,12 @@ same release version.
 
 Parameters
 ----------
-seed : {None, int}, optional
+seed : {None, int, array_like}, optional
     Random seed initializing the pseudo-random number generator.
-    Can be an integer in [0, 2**32] or ``None`` (the default).
-    If `seed` is ``None``, then ``dSFMT.RandomState`` will try to read
-    entropy from ``/dev/urandom`` (or the Windows analogue) if available to
+    Can be an integer in [0, 2**32-1], array of integers in
+    [0, 2**32-1] or ``None`` (the default). If `seed` is ``None``,
+    then ``dSFMT.RandomState`` will try to read entropy from
+    ``/dev/urandom`` (or the Windows analogue) if available to
     produce a 64-bit seed. If unavailable, the a 64-bit hash of the time
     and process ID is used.
 
@@ -136,6 +139,20 @@ segments come from the same sequence.
 >>> for i in range(10):
         rs[i].jump(i)
 
+**State and Seeding**
+The ``dsfmt.RandomState`` state vector consists of a 764 element array of
+32-bit unsigned integers plus a single integer value between 0 and 382
+indicating  the current position within the main array. The implementation
+used here augments this with a 384 element array of doubles which are used
+to efficiently access the random numbers produced by the dSFMT generator.
+
+``dsfmt.RandomState`` is seeded using either a single 32-bit unsigned integer
+or a vector of 32-bit unsigned integers.  In either case the input seed is
+used as an input (or inputs) for a hashing function, and the output of the
+hashing function is used as the initial state. Using a single 32-bit value
+for the seed can only initialize a small range of the possible initial
+state values.
+
 .. [1] Mutsuo Saito and Makoto Matsumoto, "SIMD-oriented Fast Mersenne
        Twister: a 128-bit Pseudorandom Number Generator." Monte Carlo
        and Quasi-Monte Carlo Methods 2006, Springer, pp. 607 -- 622, 2008.

randomstate/shims/mlfg-1279-861/mlfg-1279-861-shim.c

@@ -11,6 +11,11 @@ void set_seed(aug_state* state, uint64_t val)
     mlfg_seed(state->rng, val);
 }
 
+void set_seed_by_array(aug_state* state, uint64_t *vals, int count)
+{
+    mlfg_seed_by_array(state->rng, vals, count);
+}
+
 void entropy_init(aug_state* state)
 {
     uint64_t seeds[1];

randomstate/shims/mlfg-1279-861/mlfg-1279-861-shim.h

@@ -45,4 +45,6 @@ inline double random_double(aug_state* state)
 
 extern void set_seed(aug_state* state, uint64_t seed);
 
+extern void set_seed_by_array(aug_state* state, uint64_t *vals, int count);
+
 extern void entropy_init(aug_state* state);

randomstate/shims/mlfg-1279-861/mlfg-1279-861.pxi

@@ -22,6 +22,8 @@ cdef extern from "distributions.h":
 
     cdef void set_seed(aug_state* state, uint64_t seed)
 
+    cdef void set_seed_by_array(aug_state* state, uint64_t *seed, int count)
+
 ctypedef uint64_t rng_state_t
 
 ctypedef mlfg_state rng_t
@@ -76,17 +78,35 @@ same release version.
 
 Parameters
 ----------
-seed : {None, int}, optional
+seed : {None, int, array_like}, optional
     Random seed initializing the pseudo-random number generator.
-    Can be an integer in [0, 2**64] or ``None`` (the default).
-    If ``seed`` is ``None``, then ``mlfg_1279_861.RandomState`` will try to
-    read entropy from ``/dev/urandom`` (or the Windows analogue) if available
-    to produce a 64-bit seed. If unavailable, a 64-bit hash of the time
+    Can be an integer in [0, 2**64-1], array of integers in
+    [0, 2**64-1] or ``None`` (the default). If `seed` is ``None``,
+    then ``mlfg_1279_861.RandomState`` will try to read entropy from
+    ``/dev/urandom`` (or the Windows analogue) if available to
+    produce a 64-bit seed. If unavailable, a 64-bit hash of the time
     and process ID is used.
 
 Notes
 -----
 The state of the MLFG(1279,861,*) PRNG is represented by 1279 64-bit unsigned
 integers as well as a two 32-bit integers representing the location in the
 state array of the current and lagged values.
+
+**State and Seeding**
+The ``dsfmt.mlfg_1279_861`` state vector consists of a 1279 element array of
+64-bit unsigned integers plus a two integers value between 0 and 1278
+indicating  the current position and the lagged value within the main array
+required to produce the next random.
+
+``dsfmt.mlfg_1279_861`` is seeded using either a single 64-bit unsigned integer
+or a vector of 64-bit unsigned integers.  In either case the input seed is
+used as an input (or inputs) for another simple random number generator,
+Splitmix64, and the output of this PRNG function is used as the initial state.
+Using a single 64-bit value for the seed can only initialize a small range of
+the possible initial state values.  When using an array, the SplitMix64 state
+for producing the ith component of the initial state is XORd with the ith
+value of the seed array until the seed array is exhausted. When using an array
+the initial state for the SplitMix64 state is 0 so that using a single element
+array and using a single scalar value will produce the same initial state.
 """

randomstate/shims/mrg32k3a/mrg32k3a-shim.c

@@ -30,6 +30,12 @@ void set_seed(aug_state* state, uint64_t val)
     mrg32k3a_seed(state->rng, val);
 }
 
+void set_seed_by_array(aug_state* state, uint64_t *vals, int count)
+{
+    mrg32k3a_seed_by_array(state->rng, vals, count);
+}
+
+
 void init_state(aug_state* state, int64_t vals[6])
 {
     state->rng->s1[0] =  vals[0];

randomstate/shims/mrg32k3a/mrg32k3a-shim.h

@@ -39,6 +39,8 @@ inline double random_double(aug_state* state)
 
 extern void set_seed(aug_state* state, uint64_t val);
 
+extern void set_seed_by_array(aug_state* state, uint64_t *vals, int count);
+
 extern void init_state(aug_state* state, int64_t vals[6]);
 
 extern void entropy_init(aug_state* state);

randomstate/shims/mrg32k3a/mrg32k3a.pxi

@@ -23,6 +23,8 @@ cdef extern from "distributions.h":
 
     cdef void set_seed(aug_state* state, uint64_t seed)
 
+    cdef void set_seed_by_array(aug_state* state, uint64_t *seed, int count)
+
 ctypedef mrg32k3a_state rng_t
 
 ctypedef uint64_t rng_state_t
@@ -129,11 +131,12 @@ same release version.
 
 Parameters
 ----------
-seed : {None, int}, optional
+seed : {None, int, array_like}, optional
     Random seed initializing the pseudo-random number generator.
-    Can be an integer in [0, 2**64] or ``None`` (the default).
-    If `seed` is ``None``, then ``mrg32k3a.RandomState`` will try to read data
-    from ``/dev/urandom`` (or the Windows analogue) if available. If
+    Can be an integer in [0, 2**64-1], array of integers in
+    [0, 2**64-1] or ``None`` (the default). If `seed` is ``None``,
+    then ``mrg32k3a.RandomState`` will try to read data from
+    ``/dev/urandom`` (or the Windows analogue) if available. If
     unavailable, a 64-bit hash of the time and process ID is used.
 
 Notes

randomstate/shims/xorshift1024/xorshift1024-shim.c

@@ -11,6 +11,11 @@ void set_seed(aug_state* state, uint64_t seed)
     xorshift1024_seed(state->rng, seed);
 }
 
+void set_seed_by_array(aug_state* state, uint64_t *vals, int count)
+{
+    xorshift1024_seed_by_array(state->rng, vals, count);
+}
+
 void entropy_init(aug_state* state)
 {
     uint64_t seed[1];