Merge pull request #56 from bashtage/add-dtypes-uniform

ENH: Add single precision uniforms

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
 

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
 

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
 --------------

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.

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 <parallel>
+   Reading System Entropy <entropy>
+
 
 Individual Pseudo Random Number Generators
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -85,8 +97,13 @@ Individual Pseudo Random Number Generators
    PCG-64 <pcg64>
    MLFG <mlfg>
    MRG32K3A <mrg32k3a>
-   Reading System Entropy <entropy>
 
+Changes
+~~~~~~~
+.. toctree::
+   :maxdepth: 2
+
+   Change Log <change-log>
 
 Indices and tables
 ~~~~~~~~~~~~~~~~~~

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.ndarray>np.empty(size, np.float64)
         n = np.PyArray_SIZE(out_array)
         out_array_data = <double *>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 = <random_float_fill>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.ndarray>np.empty(size, np.float32)
+        n = np.PyArray_SIZE(out_array)
+        out_array_data = <float *>np.PyArray_DATA(out_array)
         with lock, nogil:
             f(state, n, out_array_data)
         return out_array

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++) {

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);
 

randomstate/interface/dSFMT/dSFMT-shim.h

@@ -39,14 +39,16 @@ 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);
 }
 
 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;

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):
         """