From e34e8e4bbd7a2f137d52579bf011b77545f8b466 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Mon, 13 Jun 2016 12:02:17 +0100 Subject: [PATCH] ENH: Add single precision uniforms Add precision choice for uniform random numbers to allow singles (float32) or doubles (float64). --- README.md | 3 + README.rst | 1 + doc/source/change-log.rst | 15 ++++ doc/source/conf.py | 4 +- doc/source/index.rst | 23 +++++- randomstate/array_utilities.pxi | 21 +++++ randomstate/distributions.c | 16 +++- randomstate/distributions.h | 8 +- randomstate/interface/dSFMT/dSFMT-shim.h | 4 +- randomstate/randomstate.pyx | 19 ++++- randomstate/tests/test_direct.py | 98 ++++++++++++++++++++---- setup.py | 2 +- 12 files changed, 187 insertions(+), 27 deletions(-) diff --git a/README.md b/README.md index 82737f44..7ed97f7d 100644 --- a/README.md +++ b/README.md @@ -53,6 +53,9 @@ version of the MT19937 generator that is especially fast at generating doubles support an additional `method` keyword argument which can be `bm` or `zig` where `bm` corresponds to the current method using the Box-Muller transformation and `zig` uses the much faster (100%+) ziggurat method. +* `random_sample` can produce either single precision (`np.float32`) or +double precision (`np.float64`, the default) using an the optional keyword +argument `dtype`. ### New Functions diff --git a/README.rst b/README.rst index 7fb7662d..bdf0a4b3 100644 --- a/README.rst +++ b/README.rst @@ -7,6 +7,7 @@ This is a library and generic interface for alternative random generators in Python and Numpy. Features +-------- - Immediate drop in replacement for NumPy's RandomState diff --git a/doc/source/change-log.rst b/doc/source/change-log.rst index a7b81707..b0f4abe7 100644 --- a/doc/source/change-log.rst +++ b/doc/source/change-log.rst @@ -2,6 +2,21 @@ Change Log ========== +Version 1.11.2 +-------------- +* Added keyword argument `dtype` to `random_sample` which allows for single + precision as well as double precision uniforms to be generated. + +.. ipython:: python + + import numpy as np + import randomstate as rs + rs.seed(23456) + rs.random_sample(3, dtype=np.float64) + + rs.seed(23456) + rs.random_sample(3, dtype=np.float32) + Version 1.11.1 -------------- diff --git a/doc/source/conf.py b/doc/source/conf.py index d70cacc1..849ecbc7 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -37,6 +37,8 @@ 'sphinx.ext.mathjax', 'sphinx.ext.napoleon', 'sphinx.ext.autosummary', + 'IPython.sphinxext.ipython_console_highlighting', + 'IPython.sphinxext.ipython_directive' ] # Add any paths that contain templates here, relative to this directory. @@ -65,7 +67,7 @@ # The short X.Y version. version = '1.11' # The full version, including alpha/beta/rc tags. -release = '1.11.1' +release = '1.11.2' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/doc/source/index.rst b/doc/source/index.rst index 43400a22..c5f34fc4 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -10,8 +10,12 @@ What's New or Different * ``randomstate.entropy.random_entropy`` provides access to the system source of randomness that is used in cryptographic applications (e.g., ``/dev/urandom`` on Unix). -* The normal generator supports a 256-step ziggurat method which is 2-6 times - faster than NumPy's ``standard_normal``. This generator can be accessed using +* The normal generator supports a 256-step Ziggurat method which is 2-6 times + faster than NumPy's ``standard_normal``. This generator can be accessed + by passing the keyword argument ``method='zig'``. +* ``random_sample`` accepts the optional keyword argument ``dtype`` which + accepts ``np.float32`` or ``np.float64`` to produce either single or + double prevision uniform random variables * For changes since the previous release, see the :ref:`change-log` .. code-block:: python @@ -69,6 +73,14 @@ generators, 'in addition' to the standard PRNG in NumPy. The included PRNGs are .. _`wiki page on Fibonacci generators`: https://en.wikipedia.org/wiki/Lagged_Fibonacci_generator .. _`MRG32K3A author's page`: http://simul.iro.umontreal.ca/ +New Features +~~~~~~~~~~~~ +.. toctree:: + :maxdepth: 2 + + Using in Parallel Applications + Reading System Entropy + Individual Pseudo Random Number Generators ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -85,8 +97,13 @@ Individual Pseudo Random Number Generators PCG-64 MLFG MRG32K3A - Reading System Entropy +Changes +~~~~~~~ +.. toctree:: + :maxdepth: 2 + + Change Log Indices and tables ~~~~~~~~~~~~~~~~~~ diff --git a/randomstate/array_utilities.pxi b/randomstate/array_utilities.pxi index 38adf8c6..e96e1625 100644 --- a/randomstate/array_utilities.pxi +++ b/randomstate/array_utilities.pxi @@ -16,8 +16,10 @@ ctypedef uint32_t (* random_uint_1_i_32)(aug_state* state, uint32_t a) nogil ctypedef int32_t (* random_int_2_i_32)(aug_state* state, int32_t a, int32_t b) nogil ctypedef int64_t (* random_int_2_i)(aug_state* state, int64_t a, int64_t b) nogil +ctypedef void (* random_float_fill)(aug_state* state, np.npy_intp count, float *out) nogil ctypedef void (* random_double_fill)(aug_state* state, np.npy_intp count, double *out) nogil + cdef np.npy_intp compute_numel(size): cdef np.npy_intp i, n = 1 if isinstance(size, tuple): @@ -613,6 +615,25 @@ cdef object double_fill(aug_state* state, void *func, object size, object lock): out_array = np.empty(size, np.float64) n = np.PyArray_SIZE(out_array) out_array_data = np.PyArray_DATA(out_array) + with lock, nogil: + f(state, n, out_array_data) + return out_array + +cdef object float_fill(aug_state* state, void *func, object size, object lock): + cdef random_float_fill f = func + cdef float out + cdef float *out_array_data + cdef np.ndarray out_array + cdef np.npy_intp n + + if size is None: + with lock: + f(state, 1, &out) + return out + else: + out_array = np.empty(size, np.float32) + n = np.PyArray_SIZE(out_array) + out_array_data = np.PyArray_DATA(out_array) with lock, nogil: f(state, n, out_array_data) return out_array \ No newline at end of file diff --git a/randomstate/distributions.c b/randomstate/distributions.c index 66d2417e..6c1a1eb5 100644 --- a/randomstate/distributions.c +++ b/randomstate/distributions.c @@ -30,13 +30,25 @@ unsigned long random_uint(aug_state* state) } -double random_standard_uniform(aug_state* state) +float random_standard_uniform32(aug_state* state) +{ + return (random_uint32(state) >> 9) * (1.0 / 8388608.0); +} + +double random_standard_uniform64(aug_state* state) { return random_double(state); } +void random_uniform_fill32(aug_state* state, npy_intp count, float *out) +{ + int i; + for (i=0; i < count; i++) { + out[i] = (random_uint32(state) >> 9) * (float)(1.0 / 8388608.0); + } +} -void random_uniform_fill(aug_state* state, npy_intp count, double *out) +void random_uniform_fill64(aug_state* state, npy_intp count, double *out) { int i; for (i=0; i < count; i++) { diff --git a/randomstate/distributions.h b/randomstate/distributions.h index e666693d..33294903 100644 --- a/randomstate/distributions.h +++ b/randomstate/distributions.h @@ -48,7 +48,9 @@ extern int64_t random_positive_int64(aug_state* state); extern int32_t random_positive_int32(aug_state* state); -extern double random_standard_uniform(aug_state* state); +extern float random_standard_uniform32(aug_state* state); + +extern double random_standard_uniform64(aug_state* state); extern double random_standard_exponential(aug_state* state); @@ -138,7 +140,9 @@ extern void random_bounded_uint8_fill(aug_state *state, uint8_t off, uint8_t rng extern void random_bounded_bool_fill(aug_state *state, npy_bool off, npy_bool rng, npy_intp cnt, npy_bool *out); -extern void random_uniform_fill(aug_state* state, npy_intp count, double *out); +extern void random_uniform_fill32(aug_state* state, npy_intp count, float *out); + +extern void random_uniform_fill64(aug_state* state, npy_intp count, double *out); extern void random_standard_exponential_fill(aug_state* state, npy_intp count, double *out); diff --git a/randomstate/interface/dSFMT/dSFMT-shim.h b/randomstate/interface/dSFMT/dSFMT-shim.h index 54d66718..719a0753 100644 --- a/randomstate/interface/dSFMT/dSFMT-shim.h +++ b/randomstate/interface/dSFMT/dSFMT-shim.h @@ -39,6 +39,7 @@ static inline double random_double_from_buffer(aug_state *state) static inline uint32_t random_uint32(aug_state* state) { + /* TODO: This can be improved to use upper bits */ double d = random_double_from_buffer(state);//dsfmt_genrand_close1_open2(state->rng); uint64_t *out = (uint64_t *)&d; return (uint32_t)(*out & 0xffffffff); @@ -46,7 +47,8 @@ static inline uint32_t random_uint32(aug_state* state) static inline uint64_t random_uint64(aug_state* state) { - double d = random_double_from_buffer(state);//dsfmt_genrand_close1_open2(state->rng); + /* TODO: This can be improved to use upper bits */ + double d = random_double_from_buffer(state); //dsfmt_genrand_close1_open2(state->rng); uint64_t out; uint64_t *tmp; tmp = (uint64_t *)&d; diff --git a/randomstate/randomstate.pyx b/randomstate/randomstate.pyx index 7c0c94aa..82da14b1 100644 --- a/randomstate/randomstate.pyx +++ b/randomstate/randomstate.pyx @@ -68,7 +68,9 @@ cdef extern from "distributions.h": cdef void entropy_init(aug_state* state) nogil - cdef double random_standard_uniform(aug_state* state) nogil + cdef float random_standard_uniform32(aug_state* state) nogil + + cdef double random_standard_uniform64(aug_state* state) nogil cdef double random_gauss(aug_state* state) nogil cdef double random_gauss_zig(aug_state* state) nogil cdef double random_gauss_zig_julia(aug_state* state) nogil @@ -114,7 +116,8 @@ cdef extern from "distributions.h": cdef void random_bounded_uint16_fill(aug_state *state, uint16_t off, uint16_t rng, intptr_t cnt, uint16_t *out) nogil cdef void random_bounded_uint8_fill(aug_state *state, uint8_t off, uint8_t rng, intptr_t cnt, uint8_t *out) nogil cdef void random_bounded_bool_fill(aug_state *state, np.npy_bool off, np.npy_bool rng, intptr_t cnt, np.npy_bool *out) nogil - cdef void random_uniform_fill(aug_state *state, intptr_t cnt, double *out) nogil + cdef void random_uniform_fill32(aug_state *state, intptr_t cnt, double *out) nogil + cdef void random_uniform_fill64(aug_state *state, intptr_t cnt, double *out) nogil cdef void random_standard_exponential_fill(aug_state* state, intptr_t count, double *out) nogil cdef void random_gauss_fill(aug_state* state, intptr_t count, double *out) nogil cdef void random_gauss_zig_julia_fill(aug_state* state, intptr_t count, double *out) nogil @@ -680,7 +683,7 @@ cdef class RandomState: self.get_state()) # Basic distributions: - def random_sample(self, size=None): + def random_sample(self, size=None, dtype=np.float64): """ random_sample(size=None) @@ -698,6 +701,9 @@ cdef class RandomState: Output shape. If the given shape is, e.g., ``(m, n, k)``, then ``m * n * k`` samples are drawn. Default is None, in which case a single value is returned. + dtype : dtype, optional + Desired dtype of the result, either np.float64 (default) + or np.float32. Returns ------- @@ -722,7 +728,12 @@ cdef class RandomState: [-1.23204345, -1.75224494]]) """ - return double_fill(&self.rng_state, &random_uniform_fill, size, self.lock) + if dtype is np.float64: + return double_fill(&self.rng_state, &random_uniform_fill64, size, self.lock) + elif dtype is np.float32: + return float_fill(&self.rng_state, &random_uniform_fill32, size, self.lock) + else: + raise ValueError('Unknown dtype') def tomaxint(self, size=None): """ diff --git a/randomstate/tests/test_direct.py b/randomstate/tests/test_direct.py index c598b7b7..013fcc95 100644 --- a/randomstate/tests/test_direct.py +++ b/randomstate/tests/test_direct.py @@ -2,6 +2,7 @@ import sys from os.path import join from unittest import TestCase +from nose import SkipTest import numpy as np from randomstate.prng.mlfg_1279_861 import mlfg_1279_861 @@ -21,15 +22,58 @@ pwd = os.path.dirname(os.path.abspath(__file__)) -def uniform_from_uint(x, bits): +def uniform32_from_uint64(x): + x = np.uint64(x) + upper = np.array(x >> np.uint64(32), dtype=np.uint32) + lower = np.uint64(0xffffffff) + lower = np.array(x & lower, dtype=np.uint32) + joined = np.column_stack([lower, upper]).ravel() + out = (joined >> np.uint32(9)) * (1.0 / 2 ** 23) + return out.astype(np.float32) + + +def uniform32_from_uint63(x): + x = np.uint64(x) + x = np.uint32(x >> np.uint64(32)) + out = (x >> np.uint32(9)) * (1.0 / 2 ** 23) + return out.astype(np.float32) + + +def uniform32_from_uint53(x): + x = np.uint64(x) + x = np.uint32(x & np.uint64(0xffffffff)) + out = (x >> np.uint32(9)) * (1.0 / 2 ** 23) + return out.astype(np.float32) + + +def uniform32_from_uint32(x): + return (x >> np.uint32(9)) * (1.0 / 2 ** 23) + + +def uniform32_from_uint(x, bits): if bits == 64: + return uniform32_from_uint64(x) + elif bits == 63: + return uniform32_from_uint63(x) + elif bits == 53: + return uniform32_from_uint53(x) + elif bits == 32: + return uniform32_from_uint32(x) + else: + raise NotImplementedError + + +def uniform_from_uint(x, bits): + if bits in (64, 63, 53): return uniform_from_uint64(x) elif bits == 32: return uniform_from_uint32(x) + def uniform_from_uint64(x): return (x >> np.uint64(11)) * (1.0 / 9007199254740992.0) + def uniform_from_uint32(x): out = np.empty(len(x) // 2) for i in range(0, len(x), 2): @@ -38,6 +82,7 @@ def uniform_from_uint32(x): out[i // 2] = (a * 67108864.0 + b) / 9007199254740992.0 return out + def uint64_from_uint63(x): out = np.empty(len(x) // 2, dtype=np.uint64) for i in range(0, len(x), 2): @@ -46,11 +91,13 @@ def uint64_from_uint63(x): out[i // 2] = a | b return out + def uniform_from_dsfmt(x): return x.view(np.double) - 1.0 + def gauss_from_uint(x, n, bits): - if bits == 64: + if bits in (64, 63): doubles = uniform_from_uint64(x) elif bits == 32: doubles = uniform_from_uint32(x) @@ -58,15 +105,16 @@ def gauss_from_uint(x, n, bits): doubles = uniform_from_dsfmt(x) gauss = [] loc = 0 + x1 = x2 = 0.0 while len(gauss) < n: r2 = 2 while r2 >= 1.0 or r2 == 0.0: - x1 = 2.0 * doubles[loc] - 1.0; - x2 = 2.0 * doubles[loc + 1] - 1.0; - r2 = x1 * x1 + x2 * x2; + x1 = 2.0 * doubles[loc] - 1.0 + x2 = 2.0 * doubles[loc + 1] - 1.0 + r2 = x1 * x1 + x2 * x2 loc += 2 - f = np.sqrt(-2.0 * np.log(r2) / r2); + f = np.sqrt(-2.0 * np.log(r2) / r2) gauss.append(f * x2) gauss.append(f * x1) @@ -74,6 +122,9 @@ def gauss_from_uint(x, n, bits): class Base(object): + dtype = np.uint64 + data2 = data1 = {} + @classmethod def setUpClass(cls): cls.RandomState = xorshift128.RandomState @@ -93,11 +144,12 @@ def _read_csv(cls, filename): def test_raw(self): rs = self.RandomState(*self.data1['seed']) - uints = rs.random_uintegers(1000, bits=self.bits) + bitsize = 64 if self.bits > 32 else self.bits + uints = rs.random_uintegers(1000, bits=bitsize) assert_equal(uints, self.data1['data']) rs = self.RandomState(*self.data2['seed']) - uints = rs.random_uintegers(1000, bits=self.bits) + uints = rs.random_uintegers(1000, bits=bitsize) assert_equal(uints, self.data2['data']) def test_double(self): @@ -105,11 +157,13 @@ def test_double(self): vals = uniform_from_uint(self.data1['data'], self.bits) uniforms = rs.random_sample(len(vals)) assert_allclose(uniforms, vals) + assert_equal(uniforms.dtype, np.float64) rs = self.RandomState(*self.data2['seed']) vals = uniform_from_uint(self.data2['data'], self.bits) uniforms = rs.random_sample(len(vals)) assert_allclose(uniforms, vals) + assert_equal(uniforms.dtype, np.float64) def test_gauss_inv(self): n = 25 @@ -123,6 +177,19 @@ def test_gauss_inv(self): assert_allclose(gauss, gauss_from_uint(self.data2['data'], n, self.bits)) + def test_32bit_uniform(self): + rs = self.RandomState(*self.data1['seed']) + vals = uniform32_from_uint(self.data1['data'], self.bits) + uniforms = rs.random_sample(len(vals), dtype=np.float32) + assert_allclose(uniforms, vals) + assert_equal(uniforms.dtype, np.float32) + + rs = self.RandomState(*self.data2['seed']) + vals = uniform32_from_uint(self.data2['data'], self.bits) + uniforms = rs.random_sample(len(vals), dtype=np.float32) + assert_allclose(uniforms, vals) + assert_equal(uniforms.dtype, np.float32) + class TestXorshift128(Base, TestCase): @classmethod @@ -132,6 +199,7 @@ def setUpClass(cls): cls.dtype = np.uint64 cls.data1 = cls._read_csv(join(pwd, './data/xorshift128-testset-1.csv')) cls.data2 = cls._read_csv(join(pwd, './data/xorshift128-testset-2.csv')) + cls.uniform32_func = uniform32_from_uint64 class TestXoroshiro128plus(Base, TestCase): @@ -143,6 +211,7 @@ def setUpClass(cls): cls.data1 = cls._read_csv(join(pwd, './data/xoroshiro128plus-testset-1.csv')) cls.data2 = cls._read_csv(join(pwd, './data/xoroshiro128plus-testset-2.csv')) + class TestXorshift1024(Base, TestCase): @classmethod def setUpClass(cls): @@ -162,6 +231,7 @@ def setUpClass(cls): cls.data1 = cls._read_csv(join(pwd, './data/randomkit-testset-1.csv')) cls.data2 = cls._read_csv(join(pwd, './data/randomkit-testset-2.csv')) + class TestPCG32(Base, TestCase): @classmethod def setUpClass(cls): @@ -171,6 +241,7 @@ def setUpClass(cls): cls.data1 = cls._read_csv(join(pwd, './data/pcg32-testset-1.csv')) cls.data2 = cls._read_csv(join(pwd, './data/pcg32-testset-2.csv')) + class TestPCG64(Base, TestCase): @classmethod def setUpClass(cls): @@ -181,7 +252,6 @@ def setUpClass(cls): cls.data2 = cls._read_csv(join(pwd, './data/pcg64-testset-2.csv')) - class TestMRG32K3A(Base, TestCase): @classmethod def setUpClass(cls): @@ -196,7 +266,7 @@ class TestMLFG(Base, TestCase): @classmethod def setUpClass(cls): cls.RandomState = mlfg_1279_861.RandomState - cls.bits = 64 + cls.bits = 63 cls.dtype = np.uint64 cls.data1 = cls._read_csv(join(pwd, './data/mlfg-testset-1.csv')) cls.data2 = cls._read_csv(join(pwd, './data/mlfg-testset-2.csv')) @@ -204,19 +274,21 @@ def setUpClass(cls): def test_raw(self): rs = self.RandomState(*self.data1['seed']) vals = uint64_from_uint63(self.data1['data']) - uints = rs.random_uintegers(len(vals), bits=self.bits) + bitsize = 64 if self.bits > 32 else self.bits + uints = rs.random_uintegers(len(vals), bits=bitsize) assert_equal(uints, vals) rs = self.RandomState(*self.data2['seed']) vals = uint64_from_uint63(self.data2['data']) - uints = rs.random_uintegers(len(vals), bits=self.bits) + uints = rs.random_uintegers(len(vals), bits=bitsize) assert_equal(uints, vals) + class TestDSFMT(Base, TestCase): @classmethod def setUpClass(cls): cls.RandomState = dsfmt.RandomState - cls.bits = 64 + cls.bits = 53 cls.dtype = np.uint64 cls.data1 = cls._read_csv(join(pwd, './data/dSFMT-testset-1.csv')) cls.data2 = cls._read_csv(join(pwd, './data/dSFMT-testset-2.csv')) diff --git a/setup.py b/setup.py index e500df9a..66e9c583 100644 --- a/setup.py +++ b/setup.py @@ -242,7 +242,7 @@ def cythonize(e, *args, **kwargs): 'Topic :: Security :: Cryptography'] setup(name='randomstate', - version='1.11.1', + version='1.11.2', classifiers=classifiers, packages=find_packages(), package_dir={'randomstate': './randomstate'},