Merge pull request #62 from bashtage/numpy-sync

ENH: Update to latest NumPy additions

Author: bashtage

randomstate/array_utilities.pxi

@@ -46,7 +46,7 @@ cdef uint64_t MAXSIZE = <uint64_t>sys.maxsize
 
 cdef int check_array_constraint(np.ndarray val, object name, constraint_type cons) except -1:
     if cons == CONS_NON_NEGATIVE:
-        if np.any(np.less(val, 0)):
+        if np.any(np.signbit(val)):
             raise ValueError(name + " < 0")
     elif cons == CONS_POSITIVE:
         if np.any(np.less_equal(val, 0)):
@@ -75,7 +75,7 @@ cdef int check_array_constraint(np.ndarray val, object name, constraint_type con
 
 cdef int check_constraint(double val, object name, constraint_type cons) except -1:
     if cons == CONS_NON_NEGATIVE:
-        if val < 0:
+        if np.signbit(val):
             raise ValueError(name + " < 0")
     elif cons == CONS_POSITIVE:
         if val <= 0:

randomstate/distributions.c

@@ -643,6 +643,11 @@ static long random_poisson_mult(aug_state *state, double lam)
     }
 }
 
+/*
+ * The transformed rejection method for generating Poisson random variables
+ * W. Hoermann
+ * Insurance: Mathematics and Economics 12, 39-45 (1993)
+ */
 #define LS2PI 0.91893853320467267
 #define TWELFTH 0.083333333333333333333333
 static long random_poisson_ptrs(aug_state *state, double lam)

randomstate/randomstate.pyx

@@ -1561,12 +1561,12 @@ cdef class RandomState:
         if method == 'bm':
             return cont(&self.rng_state, &random_normal, size, self.lock, 2,
                         loc, '', CONS_NONE,
-                        scale, 'scale', CONS_POSITIVE,
+                        scale, 'scale', CONS_NON_NEGATIVE,
                         0.0, '', CONS_NONE)
         else:
             return cont(&self.rng_state, &random_normal_zig, size, self.lock, 2,
                         loc, '', CONS_NONE,
-                        scale, 'scale', CONS_POSITIVE,
+                        scale, 'scale', CONS_NON_NEGATIVE,
                         0.0, '', CONS_NONE)
 
     def beta(self, a, b, size=None):
@@ -1659,7 +1659,7 @@ cdef class RandomState:
 
         """
         return cont(&self.rng_state, &random_exponential, size, self.lock, 1,
-                    scale, 'scale', CONS_POSITIVE,
+                    scale, 'scale', CONS_NON_NEGATIVE,
                     0.0, '', CONS_NONE,
                     0.0, '', CONS_NONE)
 
@@ -1784,12 +1784,12 @@ cdef class RandomState:
         if key == 'float64':
             return cont(&self.rng_state, &random_standard_gamma,
                         size, self.lock, 1,
-                        shape, 'shape', CONS_POSITIVE,
+                        shape, 'shape', CONS_NON_NEGATIVE,
                         0.0, '', CONS_NONE,
                         0.0, '', CONS_NONE)
         if key == 'float32':
             return cont_float(&self.rng_state, &random_standard_gamma_float,
-                              size, self.lock, shape, 'shape', CONS_POSITIVE)
+                              size, self.lock, shape, 'shape', CONS_NON_NEGATIVE)
         else:
             raise TypeError('Unsupported dtype "%s" for standard_gamma' % key)
 
@@ -1867,8 +1867,8 @@ cdef class RandomState:
 
         """
         return cont(&self.rng_state, &random_gamma, size, self.lock, 2,
-                    shape, 'shape', CONS_POSITIVE,
-                    scale, 'scale', CONS_POSITIVE,
+                    shape, 'shape', CONS_NON_NEGATIVE,
+                    scale, 'scale', CONS_NON_NEGATIVE,
                     0.0, '', CONS_NONE)
 
     def f(self, dfnum, dfden, size=None):
@@ -2616,7 +2616,7 @@ cdef class RandomState:
 
         """
         return cont(&self.rng_state, &random_weibull, size, self.lock, 1,
-                    a, 'a', CONS_POSITIVE,
+                    a, 'a', CONS_NON_NEGATIVE,
                     0.0, '', CONS_NONE,
                     0.0, '', CONS_NONE)
 
@@ -2801,7 +2801,7 @@ cdef class RandomState:
         """
         return cont(&self.rng_state, &random_laplace, size, self.lock, 2,
                     loc, 'loc', CONS_NONE,
-                    scale, 'scale', CONS_POSITIVE,
+                    scale, 'scale', CONS_NON_NEGATIVE,
                     0.0, '', CONS_NONE)
 
     def gumbel(self, loc=0.0, scale=1.0, size=None):
@@ -2918,7 +2918,7 @@ cdef class RandomState:
         """
         return cont(&self.rng_state, &random_gumbel, size, self.lock, 2,
                     loc, 'loc', CONS_NONE,
-                    scale, 'scale', CONS_POSITIVE,
+                    scale, 'scale', CONS_NON_NEGATIVE,
                     0.0, '', CONS_NONE)
 
     def logistic(self, loc=0.0, scale=1.0, size=None):
@@ -3107,7 +3107,7 @@ cdef class RandomState:
         """
         return cont(&self.rng_state, &random_lognormal, size, self.lock, 2,
                     mean, 'mean', CONS_NONE,
-                    sigma, 'sigma', CONS_POSITIVE,
+                    sigma, 'sigma', CONS_NON_NEGATIVE,
                     0.0, '', CONS_NONE)
 
     def rayleigh(self, scale=1.0, size=None):
@@ -3173,7 +3173,7 @@ cdef class RandomState:
 
         """
         return cont(&self.rng_state, &random_rayleigh, size, self.lock, 1,
-                    scale, 'scale', CONS_POSITIVE,
+                    scale, 'scale', CONS_NON_NEGATIVE,
                     0.0, '', CONS_NONE,
                     0.0, '', CONS_NONE)
 
@@ -4354,7 +4354,11 @@ cdef class RandomState:
             x_ptr = <char*><size_t>x.ctypes.data
             stride = x.strides[0]
             itemsize = x.dtype.itemsize
-            buf = np.empty_like(x[0])  # GC'd at function exit
+            # As the array x could contain python objects we use a buffer
+            # of bytes for the swaps to avoid leaving one of the objects
+            # within the buffer and erroneously decrementing it's refcount
+            # when the function exits.
+            buf = np.empty(itemsize, dtype=np.int8) # GC'd at function exit
             buf_ptr = <char*><size_t>buf.ctypes.data
             with self.lock:
                 # We trick gcc into providing a specialized implementation for

randomstate/tests/test_numpy_mt19937_regressions.py

@@ -113,5 +113,35 @@ def test_choice_sum_of_probs_tolerance(self):
             assert_(c in a)
             assert_raises(ValueError, mt19937.choice, a, p=probs*0.9)
 
+
+    def test_shuffle_of_array_of_different_length_strings(self):
+        # Test that permuting an array of different length strings
+        # will not cause a segfault on garbage collection
+        # Tests gh-7710
+        mt19937.seed(1234)
+
+        a = np.array(['a', 'a' * 1000])
+
+        for _ in range(100):
+            mt19937.shuffle(a)
+
+        # Force Garbage Collection - should not segfault.
+        import gc
+        gc.collect()
+
+    def test_shuffle_of_array_of_objects(self):
+        # Test that permuting an array of objects will not cause
+        # a segfault on garbage collection.
+        # See gh-7719
+        mt19937.seed(1234)
+        a = np.array([np.arange(1), np.arange(4)])
+
+        for _ in range(1000):
+            mt19937.shuffle(a)
+
+        # Force Garbage Collection - should not segfault.
+        import gc
+        gc.collect()
+
 if __name__ == "__main__":
     run_module_suite()

randomstate/tests/tests_numpy_mt19937.py

@@ -132,6 +132,7 @@ def test_negative_binomial(self):
         # arguments without truncation.
         self.prng.negative_binomial(0.5, 0.5)
 
+
 class TestRandint(TestCase):
 
     rfunc = random.randint
@@ -157,7 +158,6 @@ def test_rng_zero_and_extremes(self):
             lbnd = 0 if dt is np.bool_ else np.iinfo(dt).min
             ubnd = 2 if dt is np.bool_ else np.iinfo(dt).max + 1
             tgt = ubnd - 1
-
             assert_equal(self.rfunc(tgt, tgt + 1, size=1000, dtype=dt), tgt)
             tgt = lbnd
             assert_equal(self.rfunc(tgt, tgt + 1, size=1000, dtype=dt), tgt)
@@ -261,22 +261,26 @@ def test_randint(self):
 
     def test_random_integers(self):
         mt19937.seed(self.seed)
-        actual = mt19937.random_integers(-99, 99, size=(3, 2))
-        desired = np.array([[31, 3],
-                            [-52, 41],
-                            [-48, -66]])
-        assert_array_equal(actual, desired)
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", DeprecationWarning)
+            actual = mt19937.random_integers(-99, 99, size=(3, 2))
+            desired = np.array([[31, 3],
+                                [-52, 41],
+                                [-48, -66]])
+            assert_array_equal(actual, desired)
 
     def test_random_integers_max_int(self):
         # Tests whether random_integers can generate the
         # maximum allowed Python int that can be converted
         # into a C long. Previous implementations of this
-        # method have thrown an OverflowError when attemping
+        # method have thrown an OverflowError when attempting
         # to generate this integer.
-        actual = mt19937.random_integers(np.iinfo('l').max,
-                                         np.iinfo('l').max)
-        desired = np.iinfo('l').max
-        assert_equal(actual, desired)
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", DeprecationWarning)
+            actual = mt19937.random_integers(np.iinfo('l').max,
+                                             np.iinfo('l').max)
+            desired = np.iinfo('l').max
+            assert_equal(actual, desired)
 
     def test_random_integers_deprecated(self):
         with warnings.catch_warnings():
@@ -482,6 +486,10 @@ def test_exponential(self):
                             [0.68717433461363442, 1.69175666993575979]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_exponential_0(self):
+        assert_equal(mt19937.exponential(scale=0), 0)
+        assert_raises(ValueError, mt19937.exponential, scale=-0.)
+
     def test_f(self):
         mt19937.seed(self.seed)
         actual = mt19937.f(12, 77, size=(3, 2))
@@ -498,6 +506,10 @@ def test_gamma(self):
                             [31.71863275789960568, 33.30143302795922011]])
         assert_array_almost_equal(actual, desired, decimal=14)
 
+    def test_gamma_0(self):
+        assert_equal(mt19937.gamma(shape=0, scale=0), 0)
+        assert_raises(ValueError, mt19937.gamma, shape=-0., scale=-0.)
+
     def test_geometric(self):
         mt19937.seed(self.seed)
         actual = mt19937.geometric(.123456789, size=(3, 2))
@@ -514,6 +526,10 @@ def test_gumbel(self):
                             [1.10651090478803416, -0.69535848626236174]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_gumbel_0(self):
+        assert_equal(mt19937.gumbel(scale=0), 0)
+        assert_raises(ValueError, mt19937.gumbel, scale=-0.)
+
     def test_hypergeometric(self):
         mt19937.seed(self.seed)
         actual = mt19937.hypergeometric(10.1, 5.5, 14, size=(3, 2))
@@ -548,6 +564,10 @@ def test_laplace(self):
                             [-0.05391065675859356, 1.74901336242837324]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_laplace_0(self):
+        assert_equal(mt19937.laplace(scale=0), 0)
+        assert_raises(ValueError, mt19937.laplace, scale=-0.)
+
     def test_logistic(self):
         mt19937.seed(self.seed)
         actual = mt19937.logistic(loc=.123456789, scale=2.0, size=(3, 2))
@@ -556,6 +576,10 @@ def test_logistic(self):
                             [-0.21682183359214885, 2.63373365386060332]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_laplace_0(self):
+        assert_(mt19937.laplace(scale=0) in [0, 1])
+        assert_raises(ValueError, mt19937.laplace, scale=-0.)
+
     def test_lognormal(self):
         mt19937.seed(self.seed)
         actual = mt19937.lognormal(mean=.123456789, sigma=2.0, size=(3, 2))
@@ -564,6 +588,10 @@ def test_lognormal(self):
                             [65.72798501792723869, 86.84341601437161273]])
         assert_array_almost_equal(actual, desired, decimal=13)
 
+    def test_lognormal_0(self):
+        assert_equal(mt19937.lognormal(sigma=0), 1)
+        assert_raises(ValueError, mt19937.lognormal, sigma=-0.)
+
     def test_logseries(self):
         mt19937.seed(self.seed)
         actual = mt19937.logseries(p=.923456789, size=(3, 2))
@@ -654,6 +682,10 @@ def test_normal(self):
                             [4.18552478636557357, 4.46410668111310471]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_normal_0(self):
+        assert_equal(mt19937.normal(scale=0), 0)
+        assert_raises(ValueError, mt19937.normal, scale=-0.)
+
     def test_pareto(self):
         mt19937.seed(self.seed)
         actual = mt19937.pareto(a=.123456789, size=(3, 2))
@@ -701,6 +733,10 @@ def test_rayleigh(self):
                             [11.06066537006854311, 17.35468505778271009]])
         assert_array_almost_equal(actual, desired, decimal=14)
 
+    def test_rayleigh_0(self):
+        assert_equal(mt19937.rayleigh(scale=0), 0)
+        assert_raises(ValueError, mt19937.rayleigh, scale=-0.)
+
     def test_standard_cauchy(self):
         mt19937.seed(self.seed)
         actual = mt19937.standard_cauchy(size=(3, 2))
@@ -725,6 +761,10 @@ def test_standard_gamma(self):
                             [7.54838614231317084, 8.012756093271868]])
         assert_array_almost_equal(actual, desired, decimal=14)
 
+    def test_standard_gamma_0(self):
+        assert_equal(mt19937.standard_gamma(shape=0), 0)
+        assert_raises(ValueError, mt19937.standard_gamma, shape=-0.)
+
     def test_standard_normal(self):
         mt19937.seed(self.seed)
         actual = mt19937.standard_normal(size=(3, 2))
@@ -800,6 +840,10 @@ def test_weibull(self):
                             [0.67057783752390987, 1.39494046635066793]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_weibull_0(self):
+        assert_equal(mt19937.weibull(a=0), 0)
+        assert_raises(ValueError, mt19937.weibull, a=-0.)
+
     def test_zipf(self):
         mt19937.seed(self.seed)
         actual = mt19937.zipf(a=1.23, size=(3, 2))
@@ -1404,7 +1448,7 @@ def gen_random(state, out):
     def test_multinomial(self):
         def gen_random(state, out):
             out[...] = state.multinomial(10, [1/6.]*6, size=10000)
-        self.check_function(gen_random, sz=(10000,6))
+        self.check_function(gen_random, sz=(10000, 6))
 
 # See Issue #4263
 class TestSingleEltArrayInput(TestCase):