Merge pull request #115 from IntelPython/add-more-tests-for-interfaces

Add more vendored NumPy tests

Author: ndgrigorian

CHANGELOG.md

@@ -14,7 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Added third-party tests from `numpy.random` which tests the `mkl_random.interfaces.numpy_random` interface [gh-103](https://github.com/IntelPython/mkl_random/pull/103)
 
 ### Changed
-* Updates to `mkl_random` implementations to better align with newer versions of `numpy.random` [gh-103](https://github.com/IntelPython/mkl_random/pull/103)
+* Updates to `mkl_random` implementations to better align with newer versions of `numpy.random` [gh-103](https://github.com/IntelPython/mkl_random/pull/103), [gh-115](https://github.com/IntelPython/mkl_random/pull/115)
 * Made conda recipe dependency on numpy configurable through `USE_NUMPY_BASE` environment variable [gh-105](https://github.com/IntelPython/mkl_random/pull/105)
 
 ### Fixed

mkl_random/mklrand.pyx

@@ -2189,17 +2189,17 @@ cdef class _MKLRandomState:
 
         if size is not None:
             if (np.prod(size) == 0):
-                return np.empty(size, dtype=np.dtype(dtype))
+                return np.empty(size, dtype=np.dtype(_dtype))
 
-        lowbnd, highbnd, randfunc = self._choose_randint_type(dtype)
+        lowbnd, highbnd, randfunc = self._choose_randint_type(_dtype)
 
         if low < lowbnd:
             raise ValueError(
-                f"low is out of bounds for {np.dtype(dtype).name}"
+                f"low is out of bounds for {np.dtype(_dtype).name}"
             )
         if high > highbnd:
             raise ValueError(
-                f"high is out of bounds for {np.dtype(dtype).name}"
+                f"high is out of bounds for {np.dtype(_dtype).name}"
             )
         if low >= high:
             raise ValueError("low >= high")
@@ -6303,29 +6303,55 @@ cdef class _MKLRandomState:
         array([100,   0])
 
         """
-        cdef cnp.npy_intp d
+        cdef cnp.npy_intp d, sz, niter
         cdef cnp.ndarray parr "arrayObject_parr", mnarr "arrayObject_mnarr"
         cdef double *pix
         cdef int *mnix
-        cdef cnp.npy_intp sz
-
-        d = len(pvals)
-        parr = <cnp.ndarray>cnp.PyArray_ContiguousFromObject(
-            pvals, cnp.NPY_DOUBLE, 1, 1
+        cdef long ni
+
+        parr = <cnp.ndarray>cnp.PyArray_FROMANY(
+            pvals,
+            cnp.NPY_DOUBLE,
+            0,
+            1,
+            cnp.NPY_ARRAY_ALIGNED | cnp.NPY_ARRAY_C_CONTIGUOUS
         )
+        if cnp.PyArray_NDIM(parr) == 0:
+            raise TypeError("pvals must be a 1-d sequence")
+        d = cnp.PyArray_SIZE(parr)
         pix = <double*>cnp.PyArray_DATA(parr)
-
-        if kahan_sum(pix, d-1) > (1.0 + 1e-12):
-            raise ValueError("sum(pvals[:-1]) > 1.0")
-
+        if (
+            not np.all(np.greater_equal(parr, 0))
+            or not np.all(np.less_equal(parr, 1))
+        ):
+            raise ValueError("pvals < 0, pvals > 1 or pvals is NaN")
+
+        if d and kahan_sum(pix, d - 1) > (1.0 + 1e-12):
+            # When floating, but not float dtype, and close, improve the error
+            # 1.0001 works for float16 and float32
+            if (isinstance(pvals, np.ndarray)
+                    and np.issubdtype(pvals.dtype, np.floating)
+                    and pvals.dtype != float
+                    and pvals.sum() < 1.0001):
+                msg = ("sum(pvals[:-1].astype(np.float64)) > 1.0. The pvals "
+                       "array is cast to 64-bit floating point prior to "
+                       "checking the sum. Precision changes when casting may "
+                       "cause problems even if the sum of the original pvals "
+                       "is valid.")
+            else:
+                msg = "sum(pvals[:-1]) > 1.0"
+            raise ValueError(msg)
         shape = _shape_from_size(size, d)
         multin = np.zeros(shape, np.int32)
-
         mnarr = <cnp.ndarray>multin
         mnix = <int*>cnp.PyArray_DATA(mnarr)
         sz = cnp.PyArray_SIZE(mnarr)
-
-        irk_multinomial_vec(self.internal_state, sz // d, mnix, n, d, pix)
+        ni = n
+        if (ni < 0):
+            raise ValueError("n < 0")
+        # numpy#20483: Avoids divide by 0
+        niter = sz // d if d else 0
+        irk_multinomial_vec(self.internal_state, niter, mnix, n, d, pix)
 
         return multin
 
@@ -6614,11 +6640,27 @@ cdef class _MKLRandomState:
 
         """
         if isinstance(x, (int, np.integer)):
-            arr = np.arange(x)
-        else:
-            arr = np.array(x)
-        self.shuffle(arr)
-        return arr
+            # keep using long as the default here (main numpy switched to intp)
+            arr = np.arange(x, dtype=np.result_type(x, np.long))
+            self.shuffle(arr)
+            return arr
+
+        arr = np.asarray(x)
+        if arr.ndim < 1:
+            raise IndexError("x must be an integer or at least 1-dimensional")
+
+        # shuffle has fast-path for 1-d
+        if arr.ndim == 1:
+            # Return a copy if same memory
+            if np.may_share_memory(arr, x):
+                arr = np.array(arr)
+            self.shuffle(arr)
+            return arr
+
+        # Shuffle index array, dtype to ensure fast path
+        idx = np.arange(arr.shape[0], dtype=np.intp)
+        self.shuffle(idx)
+        return arr[idx]
 
 
 cdef class MKLRandomState(_MKLRandomState):

mkl_random/tests/test_random.py

@@ -224,7 +224,7 @@ class RandIntData(NamedTuple):
 def randint():
     rfunc_method = rnd.randint
     integral_dtypes = [
-        np.bool_,
+        np.bool,
         np.int8,
         np.uint8,
         np.int16,
@@ -323,8 +323,9 @@ def test_randint_respect_dtype_singleton(randint):
         assert_equal(sample.dtype, np.dtype(dt))
 
     for dt in (bool, int):
-        lbnd = 0 if dt is bool else np.iinfo(np.dtype(dt)).min
-        ubnd = 2 if dt is bool else np.iinfo(np.dtype(dt)).max + 1
+        # The legacy rng uses "long" as the default integer:
+        lbnd = 0 if dt is bool else np.iinfo("long").min
+        ubnd = 2 if dt is bool else np.iinfo("long").max + 1
 
         # gh-7284: Ensure that we get Python data types
         sample = randint.rfunc(lbnd, ubnd, dtype=dt)