ENH: Add support for filling existing arrays

Add support for core distributions for filling existing arrays
to allow parallel generation using multithreading.

closes #68

Author: bashtage

README.md

@@ -44,10 +44,13 @@ w = rnd.standard_normal(10000, method='zig')
 
   **Note**: There are _no_ plans to extend the alternative precision 
   generation to all random number types.
-  
-  
-  
 
+* Support for filling existing arrays using `out` keyword argument. Currently
+  supported in (both 32- and 64-bit outputs)
+
+    * Uniforms (`random_sample`)
+    * Exponentials (`standard_exponential`)
+    * Normals (`standard_normal`)
 
 ## Included Pseudo Random Number Generators
 
@@ -98,9 +101,9 @@ the RNG._
 `mt19937` generator is identical to `numpy.random.RandomState`, and 
 will produce an identical sequence of random numbers for a given seed.   
 * Builds and passes all tests on:
-  * Linux 32/64 bit, Python 2.7, 3.4, 3.5 (should work on 2.6 and 3.3)
+  * Linux 32/64 bit, Python 2.7, 3.4, 3.5, 3.6 (probably works on 2.6 and 3.3)
   * PC-BSD (FreeBSD) 64-bit, Python 2.7
-  * OSX  64-bit, Python 2.7
+  * OSX 64-bit, Python 2.7
   * Windows 32/64 bit (only tested on Python 2.7 and 3.5, but should 
     work on 3.3/3.4)
 
@@ -124,9 +127,9 @@ need to be smoothed.
 ## Requirements
 Building requires:
 
-  * Python (2.7, 3.4, 3.5)
-  * NumPy (1.9, 1.10, 1.11)
-  * Cython (0.22, 0.23, 0.24)
+  * Python (2.7, 3.4, 3.5, 3.6)
+  * NumPy (1.9, 1.10, 1.11, 1.12)
+  * Cython (0.22, 0.23, 0.24, 0.25)
   * tempita (0.5+), if not provided by Cython
 
 Testing requires nose (1.3+).

README.rst

@@ -1,7 +1,7 @@
 randomstate
 ===========
 
-|Travis Build Status| |Appveyor Build Status|
+|Travis Build Status| |Appveyor Build Status| |PyPI version|
 
 This is a library and generic interface for alternative random
 generators in Python and NumPy.
@@ -45,6 +45,13 @@ change.
 **Note**: There are *no* plans to extend the alternative precision
 generation to all random number types.
 
+-  Support for filling existing arrays using ``out`` keyword argument.
+   Currently supported in (both 32- and 64-bit outputs)
+
+   -  Uniforms (``random_sample``)
+   -  Exponentials (``standard_exponential``)
+   -  Normals (``standard_normal``)
+
 Included Pseudo Random Number Generators
 ----------------------------------------
 
@@ -101,7 +108,8 @@ Status
    and will produce an identical sequence of random numbers for a given
    seed.
 -  Builds and passes all tests on:
--  Linux 32/64 bit, Python 2.7, 3.4, 3.5 (should work on 2.6 and 3.3)
+-  Linux 32/64 bit, Python 2.7, 3.4, 3.5, 3.6 (probably works on 2.6 and
+   3.3)
 -  PC-BSD (FreeBSD) 64-bit, Python 2.7
 -  OSX 64-bit, Python 2.7
 -  Windows 32/64 bit (only tested on Python 2.7 and 3.5, but should work
@@ -134,9 +142,9 @@ Requirements
 
 Building requires:
 
--  Python (2.7, 3.4, 3.5)
--  NumPy (1.9, 1.10, 1.11)
--  Cython (0.22, 0.23, 0.24)
+-  Python (2.7, 3.4, 3.5, 3.6)
+-  NumPy (1.9, 1.10, 1.11, 1.12)
+-  Cython (0.22, 0.23, 0.24, 0.25)
 -  tempita (0.5+), if not provided by Cython
 
 Testing requires nose (1.3+).
@@ -264,3 +272,5 @@ NumPy's mt19937.
    :target: https://travis-ci.org/bashtage/ng-numpy-randomstate
 .. |Appveyor Build Status| image:: https://ci.appveyor.com/api/projects/status/odc5c4ukhru5xicl/branch/master?svg=true
    :target: https://ci.appveyor.com/project/bashtage/ng-numpy-randomstate/branch/master
+.. |PyPI version| image:: https://badge.fury.io/py/randomstate.svg
+   :target: https://badge.fury.io/py/randomstate

doc/source/change-log.rst

@@ -2,6 +2,17 @@
 
 Change Log
 ==========
+
+Changes since 1.11.4
+--------------------
+* Add ``out`` argument which allows filling existing arrays. This feature was
+  added to facilitate multithreaded filling of arrays using parallel random
+  number generators.
+
+  * Uniforms (:meth:`~randomstate.prng.mt19937.random_sample`)
+  * Normals (:meth:`~randomstate.prng.mt19937.standard_normal`)
+  * Standard Exponentials (:meth:`~randomstate.prng.mt19937.standard_exponential`)
+
 Version 1.11.4
 --------------
 * Fix for error in Ziggurat implementation of Normal

doc/source/conf.py

@@ -16,7 +16,7 @@
 import sys
 import os
 import shlex
-import sphinx_bootstrap_theme
+import guzzle_sphinx_theme
 
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
@@ -118,22 +118,20 @@
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
-html_theme = 'bootstrap'
-html_theme_path = sphinx_bootstrap_theme.get_html_theme_path()
 
-#on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
-#
-#if not on_rtd:  # only import and set the theme if we're building docs locally
-#    import sphinx_rtd_theme
-#    html_theme = 'sphinx_rtd_theme'
-#    html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
-
-# Theme options are theme-specific and customize the look and feel of a theme
-# further.  For a list of options available for each theme, see the
-# documentation.
-html_theme_options = {'bootswatch_theme': 'yeti',
-                      # Render the next and previous page links in navbar. (Default: true)
-                      'navbar_sidebarrel': False,}
+# Adds an HTML table visitor to apply Bootstrap table classes
+html_translator_class = 'guzzle_sphinx_theme.HTMLTranslator'
+html_theme_path = guzzle_sphinx_theme.html_theme_path()
+html_theme = 'guzzle_sphinx_theme'
+
+# Register the theme as an extension to generate a sitemap.xml
+extensions.append("guzzle_sphinx_theme")
+
+# Guzzle theme options (see theme.conf for more information)
+html_theme_options = {
+    # Set the name of the project to appear in the sidebar
+    "project_nav_name": "randomstate",
+}
 
 # Add any paths that contain custom themes here, relative to this directory.
 #html_theme_path = []

doc/source/index.rst

@@ -35,7 +35,24 @@ What's New or Different
   rs.seed(0)
   rs.random_sample(3, dtype=np.float32)
 
+* Optional ``out`` argument that allows existing arrays to be filled for
+  select distributions
+
+  * Uniforms (:meth:`~randomstate.prng.mt19937.random_sample`)
+  * Normals (:meth:`~randomstate.prng.mt19937.standard_normal`)
+  * Standard Exponentials (:meth:`~randomstate.prng.mt19937.standard_exponential`)
 
+  This allows mulththreading to fill large arrays in chunks using suitable
+  PRNGs in parallel.
+
+.. ipython:: python
+
+  import numpy as np
+  import randomstate as rs
+  existing = np.zeros(4)
+  rs.seed(0)
+  rs.random_sample(out=existing[:2])
+  print(existing)
 
 * For changes since the previous release, see the :ref:`change-log`
 

randomstate/array_fillers.pxi.in

@@ -9,24 +9,37 @@ ctypes = (('float64', 'double'),
 {{for  nptype, ctype in ctypes}}
 
 
-cdef object {{ctype}}_fill(aug_state* state, void *func, object size, object lock):
+cdef object {{ctype}}_fill(aug_state* state, void *func, object size, object lock, object out):
     cdef random_{{ctype}}_fill f = <random_{{ctype}}_fill>func
-    cdef {{ctype}} out
+    cdef {{ctype}} fout
     cdef {{ctype}} *out_array_data
     cdef np.ndarray out_array
     cdef np.npy_intp n
 
-    if size is None:
+    if size is None and out is None:
         with lock:
-            f(state, 1, &out)
-        return out
+            f(state, 1, &fout)
+        return fout
+
+    if out is not None:
+        out_array = <np.ndarray>out
+        if out_array.dtype != np.{{nptype}}:
+            raise TypeError('Supplied output array has the wrong type.'
+                             'Expected {{nptype}}, got {0}'.format(out_array.dtype))
+        if not (np.PyArray_CHKFLAGS(out_array, np.NPY_C_CONTIGUOUS) or
+                np.PyArray_CHKFLAGS(out_array, np.NPY_F_CONTIGUOUS)):
+            raise ValueError('Supplied output array is not C contiguous')
+        if size is not None:
+            # TODO: enable this !!! if tuple(size) != out_array.shape:
+            raise ValueError('size and shape of the supplied output are different')
     else:
         out_array = <np.ndarray>np.empty(size, np.{{nptype}})
-        n = np.PyArray_SIZE(out_array)
-        out_array_data = <{{ctype}} *>np.PyArray_DATA(out_array)
-        with lock, nogil:
-            f(state, n, out_array_data)
-        return out_array
+
+    n = np.PyArray_SIZE(out_array)
+    out_array_data = <{{ctype}} *>np.PyArray_DATA(out_array)
+    with lock, nogil:
+        f(state, n, out_array_data)
+    return out_array
 
 {{endfor}}
 

randomstate/randomstate.pyx

@@ -662,9 +662,9 @@ cdef class RandomState:
                 self.get_state())
 
     # Basic distributions:
-    def random_sample(self, size=None, dtype=np.float64):
+    def random_sample(self, size=None, dtype=np.float64, out=None):
         """
-        random_sample(size=None, dtype='d')
+        random_sample(size=None, dtype='d', out=None)
 
         Return random floats in the half-open interval [0.0, 1.0).
 
@@ -684,6 +684,10 @@ cdef class RandomState:
             Desired dtype of the result. All dtypes are determined by their
             name, either 'float64' or 'float32'. The default value is
             'float64'.
+        out : ndarray, optional
+            Alternative output array in which to place the result. If size is not None,
+            it must have the same shape as the provided size and must match the type of
+            the output values.
 
         Returns
         -------
@@ -710,9 +714,9 @@ cdef class RandomState:
         """
         key = np.dtype(dtype).name
         if key == 'float64':
-            return double_fill(&self.rng_state, &random_uniform_fill_double, size, self.lock)
+            return double_fill(&self.rng_state, &random_uniform_fill_double, size, self.lock, out)
         elif key == 'float32':
-            return float_fill(&self.rng_state, &random_uniform_fill_float, size, self.lock)
+            return float_fill(&self.rng_state, &random_uniform_fill_float, size, self.lock, out)
         else:
             raise TypeError('Unsupported dtype "%s" for random_sample' % key)
 
@@ -1403,9 +1407,10 @@ cdef class RandomState:
 
 
     # Complicated, continuous distributions:
-    def standard_normal(self, size=None, dtype=np.float64, method=__normal_method):
+    def standard_normal(self, size=None, dtype=np.float64, method=__normal_method,
+                        out=None):
         """
-        standard_normal(size=None, dtype='d', method='bm')
+        standard_normal(size=None, dtype='d', method='bm', out=None)
 
         Draw samples from a standard Normal distribution (mean=0, stdev=1).
 
@@ -1422,6 +1427,10 @@ cdef class RandomState:
         method : str, optional
             Either 'bm' or 'zig'. 'bm' uses the default Box-Muller transformations
             method.  'zig' uses the much faster Ziggurat method of Marsaglia and Tsang.
+        out : ndarray, optional
+            Alternative output array in which to place the result. If size is not None,
+            it must have the same shape as the provided size and must match the type of
+            the output values.
 
         Returns
         -------
@@ -1440,26 +1449,22 @@ cdef class RandomState:
         >>> s.shape
         (3, 4, 2)
 
-        Notes
-        -----
-        float32 normals are only available using the Box-Muller transformation
-
         """
         key = np.dtype(dtype).name
         if key == 'float64':
             if method == u'zig':
                 return double_fill(&self.rng_state, &random_gauss_zig_double_fill,
-                                   size, self.lock)
+                                   size, self.lock, out)
             else:
                 return double_fill(&self.rng_state, &random_gauss_fill,
-                                   size, self.lock)
+                                   size, self.lock, out)
         elif key == 'float32':
             if method == u'zig':
                 return float_fill(&self.rng_state, &random_gauss_zig_float_fill,
-                                   size, self.lock)
+                                   size, self.lock, out)
             else:
                 return float_fill(&self.rng_state, &random_gauss_fill_float,
-                                   size, self.lock)
+                                   size, self.lock, out)
         else:
             raise TypeError('Unsupported dtype "%s" for standard_normal' % key)
 
@@ -1663,9 +1668,9 @@ cdef class RandomState:
                     0.0, '', CONS_NONE,
                     0.0, '', CONS_NONE)
 
-    def standard_exponential(self, size=None, dtype=np.float64):
+    def standard_exponential(self, size=None, dtype=np.float64, out=None):
         """
-        standard_exponential(size=None, dtype=np.float64)
+        standard_exponential(size=None, dtype=np.float64, out=None)
 
         Draw samples from the standard exponential distribution.
 
@@ -1682,6 +1687,10 @@ cdef class RandomState:
             Desired dtype of the result. All dtypes are determined by their
             name, either 'float64' or 'float32'. The default value is
             'float64'.
+        out : ndarray, optional
+            Alternative output array in which to place the result. If size is not None,
+            it must have the same shape as the provided size and must match the type of
+            the output values.
 
         Returns
         -------
@@ -1699,11 +1708,11 @@ cdef class RandomState:
         if key == 'float64':
             return double_fill(&self.rng_state,
                                &random_standard_exponential_fill_double,
-                               size, self.lock)
+                               size, self.lock, out)
         elif key == 'float32':
             return float_fill(&self.rng_state,
                               &random_standard_exponential_fill_float,
-                              size, self.lock)
+                              size, self.lock, out)
         else:
             raise TypeError('Unsupported dtype "%s" for standard_exponential'
                             % key)