Merge pull request #3 from bashtage/randomgen-fixes

TST: Fix test tolerance

numpy/random/randomgen/bounded_integers.pxd.in

@@ -1,5 +1,3 @@
-from __future__ import absolute_import
-
 from libc.stdint cimport (uint8_t, uint16_t, uint32_t, uint64_t,
                           int8_t, int16_t, int32_t, int64_t, intptr_t)
 import numpy as np

numpy/random/randomgen/bounded_integers.pyx.in

@@ -1,6 +1,5 @@
 #!python
 #cython: wraparound=False, nonecheck=False, boundscheck=False, cdivision=True
-from __future__ import absolute_import
 
 import numpy as np
 cimport numpy as np

numpy/random/randomgen/common.pxd

@@ -1,7 +1,5 @@
 #cython: language_level=3
 
-from __future__ import absolute_import
-
 from libc.stdint cimport (uint8_t, uint16_t, uint32_t, uint64_t,
                           int8_t, int16_t, int32_t, int64_t, intptr_t,
                           uintptr_t)
@@ -100,14 +98,3 @@ cdef object discrete_broadcast_iii(void *func, void *state, object size, object
                                   np.ndarray a_arr, object a_name, constraint_type a_constraint,
                                   np.ndarray b_arr, object b_name, constraint_type b_constraint,
                                   np.ndarray c_arr, object c_name, constraint_type c_constraint)
-
-cdef inline void compute_complex(double *rv_r, double *rv_i, double loc_r,
-                                 double loc_i, double var_r, double var_i, double rho) nogil:
-    cdef double scale_c, scale_i, scale_r
-
-    scale_c = sqrt(1 - rho * rho)
-    scale_r = sqrt(var_r)
-    scale_i = sqrt(var_i)
-
-    rv_i[0] = loc_i + scale_i * (rho * rv_r[0] + scale_c * rv_i[0])
-    rv_r[0] = loc_r + scale_r * rv_r[0]

numpy/random/randomgen/dsfmt.pyx

@@ -10,7 +10,6 @@ except ImportError:
 import numpy as np
 cimport numpy as np
 
-from .common import interface
 from .common cimport *
 from .distributions cimport brng_t
 from .entropy import random_entropy

numpy/random/randomgen/generator.pyx

@@ -4,8 +4,7 @@ import operator
 import warnings
 
 from cpython.pycapsule cimport PyCapsule_IsValid, PyCapsule_GetPointer
-from cpython cimport (Py_INCREF, PyComplex_RealAsDouble,
-    PyComplex_ImagAsDouble, PyComplex_FromDoubles, PyFloat_AsDouble)
+from cpython cimport (Py_INCREF, PyFloat_AsDouble)
 from libc cimport string
 from libc.stdlib cimport malloc, free
 cimport numpy as np
@@ -386,10 +385,9 @@ cdef class RandomGenerator:
         """
         tomaxint(size=None)
 
-        Random integers between 0 and ``sys.maxint``, inclusive.
-
         Return a sample of uniformly distributed random integers in the interval
-        [0, ``sys.maxint``].
+        [0, ``np.iinfo(np.int).max``]. The np.int type translates to the C long
+        integer type and its precision is platform dependent.
 
         Parameters
         ----------
@@ -411,18 +409,15 @@ cdef class RandomGenerator:
 
         Examples
         --------
-        Need to construct a RandomGenerator object
-
-        >>> rg = np.random.randomgen.RandomGenerator()
+        >>> rg = np.random.randomgen.RandomGenerator() # need a RandomGenerator object
         >>> rg.tomaxint((2,2,2))
         array([[[1170048599, 1600360186], # random
                 [ 739731006, 1947757578]],
                [[1871712945,  752307660],
                 [1601631370, 1479324245]]])
-        >>> import sys
-        >>> sys.maxint
+        >>> np.iinfo(np.int).max
         2147483647
-        >>> rg.tomaxint((2,2,2)) < sys.maxint
+        >>> rg.tomaxint((2,2,2)) < np.iinfo(np.int).max
         array([[[ True,  True],
                 [ True,  True]],
                [[ True,  True],
@@ -1006,9 +1001,11 @@ cdef class RandomGenerator:
 
         Random integers of type np.int between `low` and `high`, inclusive.
 
-        Return random integers of type np.int64 from the "discrete uniform"
+        Return random integers of type np.int from the "discrete uniform"
         distribution in the closed interval [`low`, `high`].  If `high` is
-        None (the default), then results are from [1, `low`].
+        None (the default), then results are from [1, `low`]. The np.int
+        type translates to the C long integer type and its precision
+        is platform dependent.
 
         This function has been deprecated. Use randint instead.
 
@@ -1053,7 +1050,7 @@ cdef class RandomGenerator:
         4 # random
         >>> type(np.random.random_integers(5))
         <class 'numpy.int64'>
-        >>> np.random.random_integers(5, size=(3, 2))
+        >>> np.random.random_integers(5, size=(3,2))
         array([[5, 4], # random
                [3, 3],
                [4, 5]])
@@ -1267,168 +1264,6 @@ cdef class RandomGenerator:
                     0.0, '', CONS_NONE,
                     None)
 
-    def complex_normal(self, loc=0.0, gamma=1.0, relation=0.0, size=None):
-        """
-        complex_normal(loc=0.0, gamma=1.0, relation=0.0, size=None)
-
-        Draw random samples from a complex normal (Gaussian) distribution.
-
-        Parameters
-        ----------
-        loc : complex or array_like of complex
-            Mean of the distribution.
-        gamma : float, complex or array_like of float or complex
-            Variance of the distribution
-        relation : float, complex or array_like of float or complex
-            Relation between the two component normals
-        size : int or tuple of ints, optional
-            Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
-            ``m * n * k`` samples are drawn.  If size is ``None`` (default),
-            a single value is returned if ``loc``, ``gamma`` and ``relation``
-            are all scalars. Otherwise,
-            ``np.broadcast(loc, gamma, relation).size`` samples are drawn.
-
-        Returns
-        -------
-        out : ndarray or scalar
-            Drawn samples from the parameterized complex normal distribution.
-
-        See Also
-        --------
-        numpy.random.normal : random values from a real-valued normal
-            distribution
-
-        Notes
-        -----
-        **EXPERIMENTAL** Not part of official NumPy RandomState, may change until
-        formal release on PyPi.
-
-        Complex normals are generated from a bivariate normal where the
-        variance of the real component is 0.5 Re(gamma + relation), the
-        variance of the imaginary component is 0.5 Re(gamma - relation), and
-        the covariance between the two is 0.5 Im(relation).  The implied
-        covariance matrix must be positive semi-definite and so both variances
-        must be zero and the covariance must be weakly smaller than the
-        product of the two standard deviations.
-
-        References
-        ----------
-        .. [1] Wikipedia, "Complex normal distribution",
-               https://en.wikipedia.org/wiki/Complex_normal_distribution
-        .. [2] Leigh J. Halliwell, "Complex Random Variables" in "Casualty
-               Actuarial Society E-Forum", Fall 2015.
-
-        Examples
-        --------
-        Draw samples from the distribution:
-
-        >>> s = np.random.complex_normal(size=1000)
-
-        """
-        cdef np.ndarray ogamma, orelation, oloc, randoms, v_real, v_imag, rho
-        cdef double *randoms_data
-        cdef double fgamma_r, fgamma_i, frelation_r, frelation_i, frho, fvar_r, fvar_i, \
-            floc_r, floc_i, f_real, f_imag, i_r_scale, r_scale, i_scale, f_rho
-        cdef np.npy_intp i, j, n, n2
-        cdef np.broadcast it
-
-        oloc = <np.ndarray>np.PyArray_FROM_OTF(loc, np.NPY_COMPLEX128, np.NPY_ALIGNED)
-        ogamma = <np.ndarray>np.PyArray_FROM_OTF(gamma, np.NPY_COMPLEX128, np.NPY_ALIGNED)
-        orelation = <np.ndarray>np.PyArray_FROM_OTF(relation, np.NPY_COMPLEX128, np.NPY_ALIGNED)
-
-        if np.PyArray_NDIM(ogamma) == np.PyArray_NDIM(orelation) == np.PyArray_NDIM(oloc) == 0:
-            floc_r = PyComplex_RealAsDouble(loc)
-            floc_i = PyComplex_ImagAsDouble(loc)
-            fgamma_r = PyComplex_RealAsDouble(gamma)
-            fgamma_i = PyComplex_ImagAsDouble(gamma)
-            frelation_r = PyComplex_RealAsDouble(relation)
-            frelation_i = 0.5 * PyComplex_ImagAsDouble(relation)
-
-            fvar_r = 0.5 * (fgamma_r + frelation_r)
-            fvar_i = 0.5 * (fgamma_r - frelation_r)
-            if fgamma_i != 0:
-                raise ValueError('Im(gamma) != 0')
-            if fvar_i < 0:
-                raise ValueError('Re(gamma - relation) < 0')
-            if fvar_r < 0:
-                raise ValueError('Re(gamma + relation) < 0')
-            f_rho = 0.0
-            if fvar_i > 0 and fvar_r > 0:
-                f_rho = frelation_i / sqrt(fvar_i * fvar_r)
-            if f_rho > 1.0 or f_rho < -1.0:
-                raise ValueError('Im(relation) ** 2 > Re(gamma ** 2 - relation** 2)')
-
-            if size is None:
-                f_real = random_gauss_zig(self._brng)
-                f_imag = random_gauss_zig(self._brng)
-
-                compute_complex(&f_real, &f_imag, floc_r, floc_i, fvar_r, fvar_i, f_rho)
-                return PyComplex_FromDoubles(f_real, f_imag)
-
-            randoms = <np.ndarray>np.empty(size, np.complex128)
-            randoms_data = <double *>np.PyArray_DATA(randoms)
-            n = np.PyArray_SIZE(randoms)
-
-            i_r_scale = sqrt(1 - f_rho * f_rho)
-            r_scale = sqrt(fvar_r)
-            i_scale = sqrt(fvar_i)
-            j = 0
-            with self.lock, nogil:
-                for i in range(n):
-                    f_real = random_gauss_zig(self._brng)
-                    f_imag = random_gauss_zig(self._brng)
-                    randoms_data[j+1] = floc_i + i_scale * (f_rho * f_real + i_r_scale * f_imag)
-                    randoms_data[j] = floc_r + r_scale * f_real
-                    j += 2
-
-            return randoms
-
-        gpc = ogamma + orelation
-        gmc = ogamma - orelation
-        v_real = <np.ndarray>(0.5 * np.real(gpc))
-        if np.any(np.less(v_real, 0)):
-            raise ValueError('Re(gamma + relation) < 0')
-        v_imag = <np.ndarray>(0.5 * np.real(gmc))
-        if np.any(np.less(v_imag, 0)):
-            raise ValueError('Re(gamma - relation) < 0')
-        if np.any(np.not_equal(np.imag(ogamma), 0)):
-            raise ValueError('Im(gamma) != 0')
-
-        cov = 0.5 * np.imag(orelation)
-        rho = np.zeros_like(cov)
-        idx = (v_real.flat > 0) & (v_imag.flat > 0)
-        rho.flat[idx] = cov.flat[idx] / np.sqrt(v_real.flat[idx] * v_imag.flat[idx])
-        if np.any(cov.flat[~idx] != 0) or np.any(np.abs(rho) > 1):
-            raise ValueError('Im(relation) ** 2 > Re(gamma ** 2 - relation ** 2)')
-
-        if size is not None:
-            randoms = <np.ndarray>np.empty(size, np.complex128)
-        else:
-            it = np.PyArray_MultiIterNew4(oloc, v_real, v_imag, rho)
-            randoms = <np.ndarray>np.empty(it.shape, np.complex128)
-
-        randoms_data = <double *>np.PyArray_DATA(randoms)
-        n = np.PyArray_SIZE(randoms)
-
-        it = np.PyArray_MultiIterNew5(randoms, oloc, v_real, v_imag, rho)
-        with self.lock, nogil:
-            n2 = 2 * n  # Avoid compiler noise for cast
-            for i in range(n2):
-                randoms_data[i] = random_gauss_zig(self._brng)
-        with nogil:
-            j = 0
-            for i in range(n):
-                floc_r= (<double*>np.PyArray_MultiIter_DATA(it, 1))[0]
-                floc_i= (<double*>np.PyArray_MultiIter_DATA(it, 1))[1]
-                fvar_r = (<double*>np.PyArray_MultiIter_DATA(it, 2))[0]
-                fvar_i = (<double*>np.PyArray_MultiIter_DATA(it, 3))[0]
-                f_rho = (<double*>np.PyArray_MultiIter_DATA(it, 4))[0]
-                compute_complex(&randoms_data[j], &randoms_data[j+1], floc_r, floc_i, fvar_r, fvar_i, f_rho)
-                j += 2
-                np.PyArray_MultiIter_NEXT(it)
-
-        return randoms
-
     def standard_gamma(self, shape, size=None, dtype=np.float64, out=None):
         """
         standard_gamma(shape, size=None, dtype='d', out=None)
@@ -1838,7 +1673,7 @@ cdef class RandomGenerator:
 
         Draw samples from a noncentral chi-square distribution.
 
-        The noncentral :math:`\\chi^2` distribution is a generalisation of
+        The noncentral :math:`\\chi^2` distribution is a generalization of
         the :math:`\\chi^2` distribution.
 
         Parameters
@@ -2723,7 +2558,7 @@ cdef class RandomGenerator:
         >>> def logist(x, loc, scale):
         ...     return np.exp((loc-x)/scale)/(scale*(1+np.exp((loc-x)/scale))**2)
         >>> plt.plot(bins, logist(bins, loc, scale)*count.max()/\
-        ... logist(bins, loc, scale).max())
+        ...          logist(bins, loc, scale).max())
         >>> plt.show()
 
         """
@@ -2821,7 +2656,7 @@ cdef class RandomGenerator:
         >>> # values, drawn from a normal distribution.
         >>> b = []
         >>> for i in range(1000):
-        ...    a = 10. + np.random.randn(100)
+        ...    a = 10. + np.random.standard_normal(100)
         ...    b.append(np.product(a))
 
         >>> b = np.array(b) / np.min(b) # scale values to be positive
@@ -3161,6 +2996,7 @@ cdef class RandomGenerator:
 
         >>> sum(np.random.binomial(9, 0.1, 20000) == 0)/20000.
         # answer = 0.38885, or 38%.
+
         """
 
         # Uses a custom implementation since self._binomial is required
@@ -3285,7 +3121,7 @@ cdef class RandomGenerator:
         for each successive well, that is what is the probability of a
         single success after drilling 5 wells, after 6 wells, etc.?
 
-        >>> s = np.random.negative_binomial(1, 0.9, 100000)
+        >>> s = np.random.negative_binomial(1, 0.1, 100000)
         >>> for i in range(1, 11): # doctest: +SKIP
         ...    probability = sum(s<i) / 100000.
         ...    print(i, "wells drilled, probability of one success =", probability)
@@ -3700,7 +3536,7 @@ cdef class RandomGenerator:
     def multivariate_normal(self, mean, cov, size=None, check_valid='warn',
                             tol=1e-8):
         """
-        multivariate_normal(self, mean, cov, size=None, check_valid='warn', tol=1e-8)
+        multivariate_normal(mean, cov, size=None, check_valid='warn', tol=1e-8)
 
         Draw random samples from a multivariate normal distribution.
 
@@ -3867,7 +3703,7 @@ cdef class RandomGenerator:
 
         Draw samples from a multinomial distribution.
 
-        The multinomial distribution is a multivariate generalisation of the
+        The multinomial distribution is a multivariate generalization of the
         binomial distribution.  Take an experiment with one of ``p``
         possible outcomes.  An example of such an experiment is throwing a dice,
         where the outcome can be 1 through 6.  Each sample drawn from the
@@ -4271,7 +4107,6 @@ binomial = _random_generator.binomial
 bytes = _random_generator.bytes
 chisquare = _random_generator.chisquare
 choice = _random_generator.choice
-complex_normal = _random_generator.complex_normal
 dirichlet = _random_generator.dirichlet
 exponential = _random_generator.exponential
 f = _random_generator.f

numpy/random/randomgen/mt19937.pyx

@@ -11,7 +11,6 @@ except ImportError:
 import numpy as np
 cimport numpy as np
 
-from .common import interface
 from .common cimport *
 from .distributions cimport brng_t
 from .entropy import random_entropy