diff --git a/.flake8 b/.flake8 index cc9aef1..ff088a9 100644 --- a/.flake8 +++ b/.flake8 @@ -6,7 +6,7 @@ extend-ignore = D100, # missing docstring in public class: D101, - # missing docstring in public method: + # missing docstring in public function: D103, # missing docstring in public package: D104, @@ -26,6 +26,7 @@ extend-ignore = per-file-ignores = mkl_random/__init__.py: F401 mkl_random/interfaces/__init__.py: F401 + mkl_random/tests/**/*.py: D102 filename = *.py, *.pyx, *.pxi, *.pxd max_line_length = 80 diff --git a/CHANGELOG.md b/CHANGELOG.md index 09e4f04..9828788 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * Added `mkl_random.interfaces` with `mkl_random.interfaces.numpy_random` interface, which aliases `mkl_random` functionality to more strictly adhere to NumPy's API (i.e., drops arguments and functions which are not part of standard NumPy) [gh-92](https://github.com/IntelPython/mkl_random/pull/92) +* 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) + +### Fixed +* Various bugfixes including a hang in `zipf` when called with `np.nan` and size-1 1D arrays being cast to scalars [gh-103](https://github.com/IntelPython/mkl_random/pull/103) + ### Removed * Dropped support for Python 3.9 [gh-81](https://github.com/IntelPython/mkl_random/pull/81) diff --git a/mkl_random/interfaces/_numpy_random.py b/mkl_random/interfaces/_numpy_random.py index 675e1d4..30fe308 100644 --- a/mkl_random/interfaces/_numpy_random.py +++ b/mkl_random/interfaces/_numpy_random.py @@ -281,6 +281,18 @@ def exponential(self, scale=1.0, size=None): """ return super().exponential(scale=scale, size=size) + def tomaxint(self, size=None): + """ + tomaxint(size=None) + + Return a sample of uniformly distributed random integers in the + interval [0, ``np.iinfo("long").max``]. + + For full documentation refer to `numpy.random.RandomState.tomaxint`. + + """ + return super().tomaxint(size=size) + def standard_exponential(self, size=None): """ standard_exponential(size=None) @@ -314,27 +326,29 @@ def gamma(self, shape, scale=1.0, size=None): """ return super().gamma(shape=shape, scale=scale, size=size) - def f(self, dfn, dfd, size=None): + def f(self, dfnum, dfden, size=None): """ - f(dfn, dfd, size=None) + f(dfnum, dfden, size=None) Draw samples from an F distribution. For full documentation refer to `numpy.random.f`. """ - return super().f(dfn=dfn, dfd=dfd, size=size) + return super().f(dfnum=dfnum, dfden=dfden, size=size) - def noncentral_f(self, dfn, dfd, nonc, size=None): + def noncentral_f(self, dfnum, dfden, nonc, size=None): """ - noncentral_f(dfn, dfd, nonc, size=None) + noncentral_f(dfnum, dfden, nonc, size=None) Draw samples from a non-central F distribution. For full documentation refer to `numpy.random.noncentral_f`. """ - return super().noncentral_f(dfn=dfn, dfd=dfd, nonc=nonc, size=size) + return super().noncentral_f( + dfnum=dfnum, dfden=dfden, nonc=nonc, size=size + ) def chisquare(self, df, size=None): """ diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx index 17b5885..fd799df 100644 --- a/mkl_random/mklrand.pyx +++ b/mkl_random/mklrand.pyx @@ -30,16 +30,14 @@ cdef extern from "Python.h": void PyMem_Free(void* buf) double PyFloat_AsDouble(object ob) - long PyLong_AsLong(object ob) - - int PyErr_Occurred() - void PyErr_Clear() cdef extern from "numpy/npy_no_deprecated_api.h": pass cimport cpython.tuple +cimport cython cimport numpy as cnp +from libc.stdint cimport int64_t from libc.string cimport memcpy, memset @@ -482,6 +480,7 @@ if (r < 0): import operator import warnings +from collections.abc import Sequence import numpy as np @@ -1543,6 +1542,7 @@ def _brng_id_to_name(int brng_id): cdef class _MKLRandomState: cdef irk_state *internal_state cdef object lock + _poisson_lam_max = np.iinfo("l").max - np.sqrt(np.iinfo("l").max)*10 def __init__(self, seed=None, brng="MT19937"): self.internal_state = PyMem_Malloc(sizeof(irk_state)) @@ -1569,28 +1569,32 @@ cdef class _MKLRandomState: brng_token = irk_get_brng_and_stream_mkl( self.internal_state, &stream_id ) - try: - if seed is None: - with self.lock: + with self.lock: + try: + if seed is None: _errcode = irk_randomseed_mkl( self.internal_state, brng_token, stream_id ) - else: - idx = operator.index(seed) - if idx > int(2**32 - 1) or idx < 0: - raise ValueError("Seed must be between 0 and 4294967295") - with self.lock: + else: + idx = operator.index(seed) + if idx > int(2**32 - 1) or idx < 0: + raise ValueError( + "Seed must be between 0 and 4294967295" + ) irk_seed_mkl( self.internal_state, idx, brng_token, stream_id ) - except TypeError: - obj = np.asarray(seed) - if obj.dtype is not np.dtype("uint64"): - obj = obj.astype(np.int64, casting="safe") - if ((obj > int(2**32 - 1)) | (obj < 0)).any(): - raise ValueError("Seed must be between 0 and 4294967295") - obj = obj.astype("uint32", casting="unsafe") - with self.lock: + except TypeError: + obj = np.asarray(seed) + if obj.size == 0: + raise ValueError("Seed must be non-empty") + if obj.dtype is not np.dtype("uint64"): + obj = obj.astype(np.int64, casting="safe") + if obj.ndim != 1: + raise ValueError("Seed array must be 1-d") + if ((obj > int(2**32 - 1)) | (obj < 0)).any(): + raise ValueError("Seed must be between 0 and 4294967295") + obj = obj.astype("uint32", casting="unsafe", order="C") irk_seed_mkl_array( self.internal_state, cnp.PyArray_DATA(obj), @@ -1741,12 +1745,21 @@ cdef class _MKLRandomState: cdef int brng_id cdef cnp.ndarray obj "arrayObject_obj" - if isinstance(state, (tuple, list)): - state_len = len(state) - if (state_len != 2): + if isinstance(state, dict): + try: + state = (state["bit_generator"], state["state"]["mkl_stream"]) + except KeyError: + raise ValueError("state dictionary is not valid") + algorithm_name = state[0] + else: + if not isinstance(state, (tuple, list)): + raise TypeError("state must be a dict or a tuple") + with cython.boundscheck(True): + algorithm_name = state[0] + state_len = len(state) if (state_len == 3 or state_len == 5): - algo_name, key, _pos = state[:3] - if algo_name != "MT19937": + key, _pos = state[1:3] + if algorithm_name != "MT19937": raise ValueError( "The legacy state input algorithm must be 'MT19937'" ) @@ -1765,20 +1778,9 @@ cdef class _MKLRandomState: 1, 1 ) - self.seed(obj, brng = algo_name) + self.seed(obj, brng = algorithm_name) return - raise ValueError( - "The argument to set_state must be a list of 2 elements" - ) - elif isinstance(state, dict): - try: - state = (state["bit_generator"], state["state"]["mkl_stream"]) - except KeyError: - raise ValueError("state dictionary is not valid") - else: - raise TypeError("state must be a dict or a tuple.") - algorithm_name = state[0] if algorithm_name not in _brng_dict.keys(): raise ValueError( "basic number generator algorithm must be one of ['" @@ -1881,7 +1883,7 @@ cdef class _MKLRandomState: key = np.dtype(dtype).name if key not in _randint_type: - raise TypeError('Unsupported dtype "%s" for randint' % key) + raise TypeError(f'Unsupported dtype "{key}" for randint') return _randint_type[key] # generates typed random integer in [low, high] @@ -2175,6 +2177,20 @@ cdef class _MKLRandomState: high = low low = 0 + _dtype = np.dtype(dtype) if dtype is not int else np.dtype("long") + + if not _dtype.isnative: + warnings.warn("Providing a dtype with a non-native byteorder is " + "not supported. If you require platform-independent " + "byteorder, call byteswap when required.\nIn future " + "version, providing byteorder will raise a " + "ValueError", DeprecationWarning) + _dtype = _dtype.newbyteorder() + + if size is not None: + if (np.prod(size) == 0): + return np.empty(size, dtype=np.dtype(dtype)) + lowbnd, highbnd, randfunc = self._choose_randint_type(dtype) if low < lowbnd: @@ -2191,9 +2207,8 @@ cdef class _MKLRandomState: with self.lock: ret = randfunc(low, high - 1, size) - if size is None: - if dtype in (bool, int): - return dtype(ret) + if size is None and dtype in (bool, int): + return dtype(ret) return ret @@ -2311,15 +2326,19 @@ cdef class _MKLRandomState: # __index__ must return an integer by python rules. pop_size = operator.index(a.item()) except TypeError: - raise ValueError("a must be 1-dimensional or an integer") - if pop_size <= 0: - raise ValueError("a must be greater than 0") + raise ValueError("'a' must be 1-dimensional or an integer") + if pop_size <= 0 and np.prod(size) != 0: + raise ValueError( + "'a' must be greater than 0 unless no samples are taken" + ) elif a.ndim != 1: - raise ValueError("a must be 1-dimensional") + raise ValueError("'a' must be 1-dimensional") else: pop_size = a.shape[0] - if pop_size is 0: - raise ValueError("a must be non-empty") + if pop_size is 0 and np.prod(size) != 0: + raise ValueError( + "'a' cannot be empty unless no samples are taken" + ) if p is not None: d = len(p) @@ -2329,24 +2348,33 @@ cdef class _MKLRandomState: if np.issubdtype(p.dtype, np.floating): atol = max(atol, np.sqrt(np.finfo(p.dtype).eps)) - p = cnp.PyArray_ContiguousFromObject( - p, cnp.NPY_DOUBLE, 1, 1 + p = cnp.PyArray_FROM_OTF( + p, + cnp.NPY_DOUBLE, + cnp.NPY_ARRAY_ALIGNED | cnp.NPY_ARRAY_C_CONTIGUOUS ) pix = cnp.PyArray_DATA(p) if p.ndim != 1: - raise ValueError("p must be 1-dimensional") + raise ValueError("'p' must be 1-dimensional") if p.size != pop_size: - raise ValueError("a and p must have same size") + raise ValueError("'a' and 'p' must have same size") + p_sum = kahan_sum(pix, d) + if np.isnan(p_sum): + raise ValueError("probabilities contain NaN") if np.logical_or.reduce(p < 0): raise ValueError("probabilities are not non-negative") - if abs(kahan_sum(pix, d) - 1.) > atol: + if abs(p_sum - 1.) > atol: raise ValueError("probabilities do not sum to 1") - shape = size - if shape is not None: + # `shape == None` means `shape == ()`, but with scalar unpacking at the + # end + is_scalar = size is None + if not is_scalar: + shape = size size = np.prod(shape, dtype=np.intp) else: + shape = () size = 1 # Actual sampling @@ -2356,22 +2384,25 @@ cdef class _MKLRandomState: cdf /= cdf[-1] uniform_samples = self.random_sample(shape) idx = cdf.searchsorted(uniform_samples, side="right") - idx = np.asarray(idx) # searchsorted returns a scalar + idx = np.asarray(idx) + # searchsorted returns a scalar + # force cast to int for LLP64 + idx = np.asarray(idx).astype(np.long, casting="unsafe") else: idx = self.randint(0, pop_size, size=shape) else: if size > pop_size: raise ValueError("Cannot take a larger sample than " "population when 'replace=False'") + elif size < 0: + raise ValueError("Negative dimensions are not allowed") if p is not None: if np.count_nonzero(p > 0) < size: - raise ValueError("Fewer non-zero entries in p than size") + raise ValueError("Fewer non-zero entries in 'p' than size") n_uniq = 0 p = p.copy() - found = np.zeros( - tuple() if shape is None else shape, dtype=np.int64 - ) + found = np.zeros(shape, dtype=np.long) flat_found = found.ravel() while n_uniq < size: x = self.rand(size - n_uniq) @@ -2388,10 +2419,9 @@ cdef class _MKLRandomState: idx = found else: idx = self.permutation(pop_size)[:size] - if shape is not None: - idx.shape = shape + idx = idx.reshape(shape) - if shape is None and isinstance(idx, np.ndarray): + if is_scalar and isinstance(idx, np.ndarray): # In most cases a scalar will have been made an array idx = idx.item(0) @@ -2399,7 +2429,7 @@ cdef class _MKLRandomState: if a.ndim == 0: return idx - if shape is not None and idx.ndim == 0: + if not is_scalar and idx.ndim == 0: # If size == () then the user requested a 0-d array as opposed to # a scalar object when size is None. However a[idx] is always a # scalar and not an array. So this makes sure the result is an @@ -2484,24 +2514,6 @@ cdef class _MKLRandomState: cdef cnp.ndarray olow, ohigh cdef double flow, fhigh - flow = PyFloat_AsDouble(low) - fhigh = PyFloat_AsDouble(high) - if not npy_isfinite(flow) or not npy_isfinite(fhigh): - raise OverflowError("Range exceeds valid bounds") - if flow >= fhigh: - raise ValueError("low >= high") - - if not PyErr_Occurred(): - return vec_cont2_array_sc( - self.internal_state, - irk_uniform_vec, - size, - flow, - fhigh, - self.lock - ) - - PyErr_Clear() olow = cnp.PyArray_FROM_OTF( low, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED ) @@ -2509,9 +2521,26 @@ cdef class _MKLRandomState: high, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED ) + if cnp.PyArray_NDIM(olow) == cnp.PyArray_NDIM(ohigh) == 0: + flow = PyFloat_AsDouble(low) + fhigh = PyFloat_AsDouble(high) + + if not npy_isfinite(flow) or not npy_isfinite(fhigh): + raise OverflowError("Range exceeds valid bounds") + if flow >= fhigh: + raise ValueError("low >= high") + + return vec_cont2_array_sc( + self.internal_state, + irk_uniform_vec, + size, + flow, + fhigh, + self.lock + ) + if not np.all(np.isfinite(olow)) or not np.all(np.isfinite(ohigh)): raise OverflowError("Range exceeds valid bounds") - if np.any(olow >= ohigh): raise ValueError("low >= high") @@ -2714,7 +2743,7 @@ cdef class _MKLRandomState: DeprecationWarning ) - return self.randint(low, high + 1, size=size, dtype="l") + return self.randint(low, int(high) + 1, size=size, dtype="l") # Complicated, continuous distributions: def standard_normal(self, size=None, method=ICDF): @@ -2865,10 +2894,18 @@ cdef class _MKLRandomState: cdef cnp.ndarray oloc, oscale cdef double floc, fscale - floc = PyFloat_AsDouble(loc) - fscale = PyFloat_AsDouble(scale) - if not PyErr_Occurred(): - if fscale <= 0: + oloc = cnp.PyArray_FROM_OTF( + loc, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + oscale = cnp.PyArray_FROM_OTF( + scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(oloc) == cnp.PyArray_NDIM(oscale) == 0: + floc = PyFloat_AsDouble(loc) + fscale = PyFloat_AsDouble(scale) + + if np.signbit(fscale) or fscale == 0: raise ValueError("scale <= 0") method = choose_method( method, @@ -2903,14 +2940,6 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - - oloc = cnp.PyArray_FROM_OTF( - loc, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - oscale = cnp.PyArray_FROM_OTF( - scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) if np.any(np.less_equal(oscale, 0)): raise ValueError("scale <= 0") method = choose_method( @@ -2982,9 +3011,17 @@ cdef class _MKLRandomState: cdef cnp.ndarray oa, ob cdef double fa, fb - fa = PyFloat_AsDouble(a) - fb = PyFloat_AsDouble(b) - if not PyErr_Occurred(): + oa = cnp.PyArray_FROM_OTF( + a, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + ob = cnp.PyArray_FROM_OTF( + b, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(oa) == cnp.PyArray_NDIM(ob) == 0: + fa = PyFloat_AsDouble(a) + fb = PyFloat_AsDouble(b) + if fa <= 0: raise ValueError("a <= 0") if fb <= 0: @@ -2993,14 +3030,6 @@ cdef class _MKLRandomState: self.internal_state, irk_beta_vec, size, fa, fb, self.lock ) - PyErr_Clear() - - oa = cnp.PyArray_FROM_OTF( - a, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - ob = cnp.PyArray_FROM_OTF( - b, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) if np.any(np.less_equal(oa, 0)): raise ValueError("a <= 0") if np.any(np.less_equal(ob, 0)): @@ -3053,9 +3082,14 @@ cdef class _MKLRandomState: cdef cnp.ndarray oscale cdef double fscale - fscale = PyFloat_AsDouble(scale) - if not PyErr_Occurred(): - if fscale <= 0: + oscale = cnp.PyArray_FROM_OTF( + scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(oscale) == 0: + fscale = PyFloat_AsDouble(scale) + + if np.signbit(fscale) or fscale == 0: raise ValueError("scale <= 0") return vec_cont1_array_sc( self.internal_state, @@ -3065,17 +3099,59 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - - oscale = cnp.PyArray_FROM_OTF( - scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - if np.any(np.less_equal(oscale, 0.0)): + if np.any(np.signbit(oscale) | (oscale == 0)): raise ValueError("scale <= 0") return vec_cont1_array( self.internal_state, irk_exponential_vec, size, oscale, self.lock ) + def tomaxint(self, size=None): + """ + tomaxint(size=None) + + Return a sample of uniformly distributed random integers in the interval + [0, ``np.iinfo("long").max``]. + + Parameters + ---------- + 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. Default is None, in which case a + single value is returned. + + Returns + ------- + out : ndarray + Drawn samples, with shape `size`. + + See Also + -------- + randint : Uniform sampling over a given half-open interval of integers. + random_integers : Uniform sampling over a given closed interval of + integers. + + Examples + -------- + >>> RS = mkl_random.MKLRandomState() # need a MKLRandomState object + >>> RS.tomaxint((2,2,2)) + array([[[1170048599, 1600360186], + [ 739731006, 1947757578]], + [[1871712945, 752307660], + [1601631370, 1479324245]]]) + >>> import sys + >>> sys.maxint + 2147483647 + >>> RS.tomaxint((2,2,2)) < sys.maxint + array([[[ True, True], + [ True, True]], + [[ True, True], + [ True, True]]], dtype=bool) + + """ + return vec_long_disc0_array( + self.internal_state, irk_long_vec, size, self.lock + ) + def standard_exponential(self, size=None): """ standard_exponential(size=None) @@ -3179,9 +3255,13 @@ cdef class _MKLRandomState: cdef cnp.ndarray oshape cdef double fshape - fshape = PyFloat_AsDouble(shape) - if not PyErr_Occurred(): - if fshape <= 0: + oshape = cnp.PyArray_FROM_OTF( + shape, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + if cnp.PyArray_NDIM(oshape) == 0: + fshape = PyFloat_AsDouble(shape) + + if np.signbit(fshape) or fshape == 0: raise ValueError("shape <= 0") return vec_cont1_array_sc( self.internal_state, @@ -3191,11 +3271,7 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - oshape = cnp.PyArray_FROM_OTF( - shape, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - if np.any(np.less_equal(oshape, 0.0)): + if np.any(np.signbit(oshape) | (oshape == 0)): raise ValueError("shape <= 0") return vec_cont1_array( self.internal_state, @@ -3279,12 +3355,20 @@ cdef class _MKLRandomState: cdef cnp.ndarray oshape, oscale cdef double fshape, fscale - fshape = PyFloat_AsDouble(shape) - fscale = PyFloat_AsDouble(scale) - if not PyErr_Occurred(): - if fshape <= 0: + oshape = cnp.PyArray_FROM_OTF( + shape, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + oscale = cnp.PyArray_FROM_OTF( + scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(oshape) == cnp.PyArray_NDIM(oscale) == 0: + fshape = PyFloat_AsDouble(shape) + fscale = PyFloat_AsDouble(scale) + + if np.signbit(fshape) or fshape == 0: raise ValueError("shape <= 0") - if fscale <= 0: + if np.signbit(fscale) or fscale == 0: raise ValueError("scale <= 0") return vec_cont2_array_sc( self.internal_state, @@ -3295,16 +3379,9 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - oshape = cnp.PyArray_FROM_OTF( - shape, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - oscale = cnp.PyArray_FROM_OTF( - scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - if np.any(np.less_equal(oshape, 0.0)): + if np.any(np.signbit(oshape) | (oshape == 0)): raise ValueError("shape <= 0") - if np.any(np.less_equal(oscale, 0.0)): + if np.any(np.signbit(oscale) | (oscale == 0)): raise ValueError("scale <= 0") return vec_cont2_array( self.internal_state, irk_gamma_vec, size, oshape, oscale, self.lock @@ -3395,9 +3472,17 @@ cdef class _MKLRandomState: cdef cnp.ndarray odfnum, odfden cdef double fdfnum, fdfden - fdfnum = PyFloat_AsDouble(dfnum) - fdfden = PyFloat_AsDouble(dfden) - if not PyErr_Occurred(): + odfnum = cnp.PyArray_FROM_OTF( + dfnum, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + odfden = cnp.PyArray_FROM_OTF( + dfden, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(odfnum) == cnp.PyArray_NDIM(odfden) == 0: + fdfnum = PyFloat_AsDouble(dfnum) + fdfden = PyFloat_AsDouble(dfden) + if fdfnum <= 0: raise ValueError("shape <= 0") if fdfden <= 0: @@ -3406,14 +3491,6 @@ cdef class _MKLRandomState: self.internal_state, irk_f_vec, size, fdfnum, fdfden, self.lock ) - PyErr_Clear() - - odfnum = cnp.PyArray_FROM_OTF( - dfnum, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - odfden = cnp.PyArray_FROM_OTF( - dfden, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) if np.any(np.less_equal(odfnum, 0.0)): raise ValueError("dfnum <= 0") if np.any(np.less_equal(odfden, 0.0)): @@ -3490,10 +3567,26 @@ cdef class _MKLRandomState: cdef cnp.ndarray odfnum, odfden, ononc cdef double fdfnum, fdfden, fnonc - fdfnum = PyFloat_AsDouble(dfnum) - fdfden = PyFloat_AsDouble(dfden) - fnonc = PyFloat_AsDouble(nonc) - if not PyErr_Occurred(): + odfnum = cnp.PyArray_FROM_OTF( + dfnum, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + odfden = cnp.PyArray_FROM_OTF( + dfden, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + ononc = cnp.PyArray_FROM_OTF( + nonc, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if ( + cnp.PyArray_NDIM(odfnum) + == cnp.PyArray_NDIM(odfden) + == cnp.PyArray_NDIM(ononc) + == 0 + ): + fdfnum = PyFloat_AsDouble(dfnum) + fdfden = PyFloat_AsDouble(dfden) + fnonc = PyFloat_AsDouble(nonc) + if fdfnum <= 1: raise ValueError("dfnum <= 1") if fdfden <= 0: @@ -3510,18 +3603,6 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - - odfnum = cnp.PyArray_FROM_OTF( - dfnum, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - odfden = cnp.PyArray_FROM_OTF( - dfden, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - ononc = cnp.PyArray_FROM_OTF( - nonc, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - if np.any(np.less_equal(odfnum, 1.0)): raise ValueError("dfnum <= 1") if np.any(np.less_equal(odfden, 0.0)): @@ -3604,19 +3685,19 @@ cdef class _MKLRandomState: cdef cnp.ndarray odf cdef double fdf - fdf = PyFloat_AsDouble(df) - if not PyErr_Occurred(): + odf = cnp.PyArray_FROM_OTF( + df, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(odf) == 0: + fdf = PyFloat_AsDouble(df) + if fdf <= 0: raise ValueError("df <= 0") return vec_cont1_array_sc( self.internal_state, irk_chisquare_vec, size, fdf, self.lock ) - PyErr_Clear() - - odf = cnp.PyArray_FROM_OTF( - df, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) if np.any(np.less_equal(odf, 0.0)): raise ValueError("df <= 0") return vec_cont1_array( @@ -3703,9 +3784,17 @@ cdef class _MKLRandomState: cdef cnp.ndarray odf, ononc cdef double fdf, fnonc - fdf = PyFloat_AsDouble(df) - fnonc = PyFloat_AsDouble(nonc) - if not PyErr_Occurred(): + odf = cnp.PyArray_FROM_OTF( + df, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + ononc = cnp.PyArray_FROM_OTF( + nonc, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(odf) == cnp.PyArray_NDIM(ononc) == 0: + fdf = PyFloat_AsDouble(df) + fnonc = PyFloat_AsDouble(nonc) + if fdf <= 0: raise ValueError("df <= 0") if fnonc < 0: @@ -3719,14 +3808,6 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - - odf = cnp.PyArray_FROM_OTF( - df, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - ononc = cnp.PyArray_FROM_OTF( - nonc, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) if np.any(np.less_equal(odf, 0.0)): raise ValueError("df <= 0") if np.any(np.less(ononc, 0.0)): @@ -3897,19 +3978,19 @@ cdef class _MKLRandomState: cdef cnp.ndarray odf cdef double fdf - fdf = PyFloat_AsDouble(df) - if not PyErr_Occurred(): + odf = cnp.PyArray_FROM_OTF( + df, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(odf) == 0: + fdf = PyFloat_AsDouble(df) + if fdf <= 0: raise ValueError("df <= 0") return vec_cont1_array_sc( self.internal_state, irk_standard_t_vec, size, fdf, self.lock ) - PyErr_Clear() - - odf = cnp.PyArray_FROM_OTF( - df, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) if np.any(np.less_equal(odf, 0.0)): raise ValueError("df <= 0") return vec_cont1_array( @@ -3996,9 +4077,17 @@ cdef class _MKLRandomState: cdef cnp.ndarray omu, okappa cdef double fmu, fkappa - fmu = PyFloat_AsDouble(mu) - fkappa = PyFloat_AsDouble(kappa) - if not PyErr_Occurred(): + omu = cnp.PyArray_FROM_OTF( + mu, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + okappa = cnp.PyArray_FROM_OTF( + kappa, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(omu) == cnp.PyArray_NDIM(okappa) == 0: + fmu = PyFloat_AsDouble(mu) + fkappa = PyFloat_AsDouble(kappa) + if fkappa < 0: raise ValueError("kappa < 0") return vec_cont2_array_sc( @@ -4006,18 +4095,10 @@ cdef class _MKLRandomState: irk_vonmises_vec, size, fmu, - kappa, + fkappa, self.lock ) - PyErr_Clear() - - omu = cnp.PyArray_FROM_OTF( - mu, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - okappa = cnp.PyArray_FROM_OTF( - kappa, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) if np.any(np.less(okappa, 0.0)): raise ValueError("kappa < 0") return vec_cont2_array( @@ -4114,19 +4195,19 @@ cdef class _MKLRandomState: cdef cnp.ndarray oa cdef double fa - fa = PyFloat_AsDouble(a) - if not PyErr_Occurred(): + oa = cnp.PyArray_FROM_OTF( + a, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(oa) == 0: + fa = PyFloat_AsDouble(a) + if fa <= 0: raise ValueError("a <= 0") return vec_cont1_array_sc( self.internal_state, irk_pareto_vec, size, fa, self.lock ) - PyErr_Clear() - - oa = cnp.PyArray_FROM_OTF( - a, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) if np.any(np.less_equal(oa, 0.0)): raise ValueError("a <= 0") return vec_cont1_array( @@ -4227,20 +4308,20 @@ cdef class _MKLRandomState: cdef cnp.ndarray oa cdef double fa - fa = PyFloat_AsDouble(a) - if not PyErr_Occurred(): - if fa <= 0: + oa = cnp.PyArray_FROM_OTF( + a, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(oa) == 0: + fa = PyFloat_AsDouble(a) + + if np.signbit(fa) or fa == 0.0: raise ValueError("a <= 0") return vec_cont1_array_sc( self.internal_state, irk_weibull_vec, size, fa, self.lock ) - PyErr_Clear() - - oa = cnp.PyArray_FROM_OTF( - a, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - if np.any(np.less_equal(oa, 0.0)): + if np.any(np.signbit(oa) | np.equal(oa, 0.0)): raise ValueError("a <= 0") return vec_cont1_array( self.internal_state, irk_weibull_vec, size, oa, self.lock @@ -4343,20 +4424,20 @@ cdef class _MKLRandomState: cdef cnp.ndarray oa cdef double fa - fa = PyFloat_AsDouble(a) - if not PyErr_Occurred(): - if fa <= 0: + oa = cnp.PyArray_FROM_OTF( + a, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(oa) == 0: + fa = PyFloat_AsDouble(a) + + if np.signbit(fa) or fa == 0.0: raise ValueError("a <= 0") return vec_cont1_array_sc( self.internal_state, irk_power_vec, size, fa, self.lock ) - PyErr_Clear() - - oa = cnp.PyArray_FROM_OTF( - a, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - if np.any(np.less_equal(oa, 0.0)): + if np.any(np.signbit(oa) | np.equal(oa, 0.0)): raise ValueError("a <= 0") return vec_cont1_array( self.internal_state, irk_power_vec, size, oa, self.lock @@ -4443,10 +4524,16 @@ cdef class _MKLRandomState: cdef cnp.ndarray oloc, oscale cdef double floc, fscale - floc = PyFloat_AsDouble(loc) - fscale = PyFloat_AsDouble(scale) - if not PyErr_Occurred(): - if fscale <= 0: + oloc = cnp.PyArray_FROM_OTF(loc, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED) + oscale = cnp.PyArray_FROM_OTF( + scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(oloc) == cnp.PyArray_NDIM(oscale) == 0: + floc = PyFloat_AsDouble(loc) + fscale = PyFloat_AsDouble(scale) + + if np.signbit(fscale) or fscale == 0.0: raise ValueError("scale <= 0") return vec_cont2_array_sc( self.internal_state, @@ -4457,12 +4544,7 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - oloc = cnp.PyArray_FROM_OTF(loc, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED) - oscale = cnp.PyArray_FROM_OTF( - scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - if np.any(np.less_equal(oscale, 0.0)): + if np.any(np.signbit(oscale) | np.equal(oscale, 0.0)): raise ValueError("scale <= 0") return vec_cont2_array( self.internal_state, irk_laplace_vec, size, oloc, oscale, self.lock @@ -4582,10 +4664,16 @@ cdef class _MKLRandomState: cdef cnp.ndarray oloc, oscale cdef double floc, fscale - floc = PyFloat_AsDouble(loc) - fscale = PyFloat_AsDouble(scale) - if not PyErr_Occurred(): - if fscale <= 0: + oloc = cnp.PyArray_FROM_OTF(loc, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED) + oscale = cnp.PyArray_FROM_OTF( + scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(oloc) == cnp.PyArray_NDIM(oscale) == 0: + floc = PyFloat_AsDouble(loc) + fscale = PyFloat_AsDouble(scale) + + if np.signbit(fscale) or fscale == 0.0: raise ValueError("scale <= 0") return vec_cont2_array_sc( self.internal_state, @@ -4596,12 +4684,7 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - oloc = cnp.PyArray_FROM_OTF(loc, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED) - oscale = cnp.PyArray_FROM_OTF( - scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - if np.any(np.less_equal(oscale, 0.0)): + if np.any(np.signbit(oscale) | np.equal(oscale, 0.0)): raise ValueError("scale <= 0") return vec_cont2_array( self.internal_state, irk_gumbel_vec, size, oloc, oscale, self.lock @@ -4682,10 +4765,16 @@ cdef class _MKLRandomState: cdef cnp.ndarray oloc, oscale cdef double floc, fscale - floc = PyFloat_AsDouble(loc) - fscale = PyFloat_AsDouble(scale) - if not PyErr_Occurred(): - if fscale <= 0: + oloc = cnp.PyArray_FROM_OTF(loc, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED) + oscale = cnp.PyArray_FROM_OTF( + scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(oloc) == cnp.PyArray_NDIM(oscale) == 0: + floc = PyFloat_AsDouble(loc) + fscale = PyFloat_AsDouble(scale) + + if np.signbit(fscale) or fscale == 0.0: raise ValueError("scale <= 0") return vec_cont2_array_sc( self.internal_state, @@ -4696,12 +4785,7 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - oloc = cnp.PyArray_FROM_OTF(loc, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED) - oscale = cnp.PyArray_FROM_OTF( - scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - if np.any(np.less_equal(oscale, 0.0)): + if np.any(np.signbit(oscale) | np.equal(oscale, 0.0)): raise ValueError("scale <= 0") return vec_cont2_array( self.internal_state, @@ -4822,11 +4906,18 @@ cdef class _MKLRandomState: cdef cnp.ndarray omean, osigma cdef double fmean, fsigma - fmean = PyFloat_AsDouble(mean) - fsigma = PyFloat_AsDouble(sigma) + omean = cnp.PyArray_FROM_OTF( + mean, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + osigma = cnp.PyArray_FROM_OTF( + sigma, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) - if not PyErr_Occurred(): - if fsigma <= 0: + if cnp.PyArray_NDIM(omean) == cnp.PyArray_NDIM(osigma) == 0: + fmean = PyFloat_AsDouble(mean) + fsigma = PyFloat_AsDouble(sigma) + + if np.signbit(fsigma) or fsigma == 0.0: raise ValueError("sigma <= 0") method = choose_method( method, [ICDF, BOXMULLER], _method_alias_dict_gaussian_short @@ -4850,15 +4941,7 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - - omean = cnp.PyArray_FROM_OTF( - mean, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - osigma = cnp.PyArray_FROM_OTF( - sigma, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - if np.any(np.less_equal(osigma, 0.0)): + if np.any(np.signbit(osigma) | np.equal(osigma, 0.0)): raise ValueError("sigma <= 0.0") method = choose_method( @@ -4943,21 +5026,20 @@ cdef class _MKLRandomState: cdef cnp.ndarray oscale cdef double fscale - fscale = PyFloat_AsDouble(scale) + oscale = cnp.PyArray_FROM_OTF( + scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(oscale) == 0: + fscale = PyFloat_AsDouble(scale) - if not PyErr_Occurred(): - if fscale <= 0: + if np.signbit(fscale) or fscale == 0.0: raise ValueError("scale <= 0") return vec_cont1_array_sc( self.internal_state, irk_rayleigh_vec, size, fscale, self.lock ) - PyErr_Clear() - - oscale = cnp.PyArray_FROM_OTF( - scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - if np.any(np.less_equal(oscale, 0.0)): + if np.any(np.signbit(oscale) | np.equal(oscale, 0.0)): raise ValueError("scale <= 0.0") return vec_cont1_array( self.internal_state, irk_rayleigh_vec, size, oscale, self.lock @@ -5028,9 +5110,17 @@ cdef class _MKLRandomState: cdef cnp.ndarray omean, oscale cdef double fmean, fscale - fmean = PyFloat_AsDouble(mean) - fscale = PyFloat_AsDouble(scale) - if not PyErr_Occurred(): + omean = cnp.PyArray_FROM_OTF( + mean, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + oscale = cnp.PyArray_FROM_OTF( + scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(omean) == cnp.PyArray_NDIM(oscale) == 0: + fmean = PyFloat_AsDouble(mean) + fscale = PyFloat_AsDouble(scale) + if fmean <= 0: raise ValueError("mean <= 0") if fscale <= 0: @@ -5044,13 +5134,6 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - omean = cnp.PyArray_FROM_OTF( - mean, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - oscale = cnp.PyArray_FROM_OTF( - scale, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) if np.any(np.less_equal(omean, 0.0)): raise ValueError("mean <= 0.0") elif np.any(np.less_equal(oscale, 0.0)): @@ -5123,10 +5206,26 @@ cdef class _MKLRandomState: cdef cnp.ndarray oleft, omode, oright cdef double fleft, fmode, fright - fleft = PyFloat_AsDouble(left) - fright = PyFloat_AsDouble(right) - fmode = PyFloat_AsDouble(mode) - if not PyErr_Occurred(): + oleft = cnp.PyArray_FROM_OTF( + left, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + omode = cnp.PyArray_FROM_OTF( + mode, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + oright = cnp.PyArray_FROM_OTF( + right, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if ( + cnp.PyArray_NDIM(oleft) + == cnp.PyArray_NDIM(omode) + == cnp.PyArray_NDIM(oright) + == 0 + ): + fleft = PyFloat_AsDouble(left) + fright = PyFloat_AsDouble(right) + fmode = PyFloat_AsDouble(mode) + if fleft > fmode: raise ValueError("left > mode") if fmode > fright: @@ -5143,17 +5242,6 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - oleft = cnp.PyArray_FROM_OTF( - left, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - omode = cnp.PyArray_FROM_OTF( - mode, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - oright = cnp.PyArray_FROM_OTF( - right, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) - if np.any(np.greater(oleft, omode)): raise ValueError("left > mode") if np.any(np.greater(omode, oright)): @@ -5254,12 +5342,20 @@ cdef class _MKLRandomState: """ cdef cnp.ndarray on, op - cdef long ln + cdef int64_t ln cdef double fp - fp = PyFloat_AsDouble(p) - ln = PyLong_AsLong(n) - if not PyErr_Occurred(): + on = cnp.PyArray_FROM_OTF( + n, cnp.NPY_INTP, cnp.NPY_ARRAY_IN_ARRAY + ) + op = cnp.PyArray_FROM_OTF( + p, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_IN_ARRAY + ) + + if cnp.PyArray_NDIM(on) == cnp.PyArray_NDIM(op) == 0: + fp = PyFloat_AsDouble(p) + ln = n + if ln < 0: raise ValueError("n < 0") if fp < 0: @@ -5268,33 +5364,24 @@ cdef class _MKLRandomState: raise ValueError("p > 1") elif np.isnan(fp): raise ValueError("p is nan") - if n > int(2**31-1): + if ln > int(2**31-1): raise ValueError("n > 2147483647") - else: - return vec_discnp_array_sc( - self.internal_state, - irk_binomial_vec, - size, - ln, - fp, - self.lock - ) - - PyErr_Clear() + return vec_discnp_array_sc( + self.internal_state, + irk_binomial_vec, + size, + ln, + fp, + self.lock + ) - on = cnp.PyArray_FROM_OTF( - n, cnp.NPY_LONG, cnp.NPY_ARRAY_IN_ARRAY - ) - op = cnp.PyArray_FROM_OTF( - p, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_IN_ARRAY - ) if np.any(np.less(n, 0)): raise ValueError("n < 0") if np.any(np.less(op, 0)): raise ValueError("p < 0") if np.any(np.greater(op, 1)): raise ValueError("p > 1") - if np.any(np.greater(n, int(2**31-1))): + if np.any(np.greater(n, int(2**31 - 1))): raise ValueError("n > 2147483647") on = on.astype(np.int32, casting="unsafe") @@ -5377,9 +5464,17 @@ cdef class _MKLRandomState: cdef double fn cdef double fp - fp = PyFloat_AsDouble(p) - fn = PyFloat_AsDouble(n) - if not PyErr_Occurred(): + on = cnp.PyArray_FROM_OTF( + n, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_IN_ARRAY + ) + op = cnp.PyArray_FROM_OTF( + p, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_IN_ARRAY + ) + + if cnp.PyArray_NDIM(op) == cnp.PyArray_NDIM(on) == 0: + fp = PyFloat_AsDouble(p) + fn = PyFloat_AsDouble(n) + if fn <= 0: raise ValueError("n <= 0") if fp < 0: @@ -5395,14 +5490,6 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - - on = cnp.PyArray_FROM_OTF( - n, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_IN_ARRAY - ) - op = cnp.PyArray_FROM_OTF( - p, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_IN_ARRAY - ) if np.any(np.less_equal(n, 0)): raise ValueError("n <= 0") if np.any(np.less(p, 0)): @@ -5483,13 +5570,17 @@ cdef class _MKLRandomState: """ cdef cnp.ndarray olam cdef double flam - poisson_lam_max = np.iinfo("l").max - np.sqrt(np.iinfo("l").max)*10 - flam = PyFloat_AsDouble(lam) - if not PyErr_Occurred(): + olam = cnp.PyArray_FROM_OTF( + lam, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_IN_ARRAY + ) + + if cnp.PyArray_NDIM(olam) == 0: + flam = PyFloat_AsDouble(lam) + if lam < 0: raise ValueError("lam < 0") - if lam > poisson_lam_max: + if lam > self._poisson_lam_max: raise ValueError("lam value too large") method = choose_method( method, [POISNORM, PTPE], _method_alias_dict_poisson @@ -5511,14 +5602,9 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - - olam = cnp.PyArray_FROM_OTF( - lam, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_IN_ARRAY - ) if np.any(np.less(olam, 0)): raise ValueError("lam < 0") - if np.any(np.greater(olam, poisson_lam_max)): + if np.any(np.greater(olam, self._poisson_lam_max)): raise ValueError("lam value too large.") method = choose_method( method, [POISNORM, PTPE], _method_alias_dict_poisson @@ -5615,21 +5701,22 @@ cdef class _MKLRandomState: cdef cnp.ndarray oa cdef double fa - fa = PyFloat_AsDouble(a) - if not PyErr_Occurred(): - if fa <= 1.0: - raise ValueError("a <= 1.0") + oa = cnp.PyArray_FROM_OTF( + a, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_IN_ARRAY + ) + + if cnp.PyArray_NDIM(oa) == 0: + fa = PyFloat_AsDouble(a) + + if not fa > 1.0: + raise ValueError("a <= 1.0 or a is NaN") return vec_long_discd_array_sc( self.internal_state, irk_zipf_long_vec, size, fa, self.lock ) - PyErr_Clear() - - oa = cnp.PyArray_FROM_OTF( - a, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_IN_ARRAY - ) - if np.any(np.less_equal(oa, 1.0)): - raise ValueError("a <= 1.0") + # NaN is not greater than 1.0 + if not np.all(np.greater(oa, 1.0)): + raise ValueError("a <= 1.0 or a contains NaNs") return vec_long_discd_array( self.internal_state, irk_zipf_long_vec, size, oa, self.lock ) @@ -5683,8 +5770,13 @@ cdef class _MKLRandomState: cdef cnp.ndarray op cdef double fp - fp = PyFloat_AsDouble(p) - if not PyErr_Occurred(): + op = cnp.PyArray_FROM_OTF( + p, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_IN_ARRAY + ) + + if cnp.PyArray_NDIM(op) == 0: + fp = PyFloat_AsDouble(p) + if fp <= 0.0: raise ValueError("p <= 0.0") if fp > 1.0: @@ -5693,13 +5785,8 @@ cdef class _MKLRandomState: self.internal_state, irk_geometric_vec, size, fp, self.lock ) - PyErr_Clear() - - op = cnp.PyArray_FROM_OTF( - p, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_IN_ARRAY - ) if np.any(np.less_equal(op, 0.0)): - raise ValueError("p < 0.0") + raise ValueError("p <= 0.0") if np.any(np.greater(op, 1.0)): raise ValueError("p > 1.0") return vec_discd_array( @@ -5793,12 +5880,28 @@ cdef class _MKLRandomState: """ cdef cnp.ndarray ongood, onbad, onsample, otot - cdef long lngood, lnbad, lnsample, lntot + cdef int64_t lngood, lnbad, lnsample, lntot + + ongood = cnp.PyArray_FROM_OTF( + ngood, cnp.NPY_INT64, cnp.NPY_ARRAY_ALIGNED + ) + onbad = cnp.PyArray_FROM_OTF( + nbad, cnp.NPY_INT64, cnp.NPY_ARRAY_ALIGNED + ) + onsample = cnp.PyArray_FROM_OTF( + nsample, cnp.NPY_INT64, cnp.NPY_ARRAY_ALIGNED + ) + + if ( + cnp.PyArray_NDIM(ongood) + == cnp.PyArray_NDIM(onbad) + == cnp.PyArray_NDIM(onsample) + == 0 + ): + lngood = ngood + lnbad = nbad + lnsample = nsample - lngood = PyLong_AsLong(ngood) - lnbad = PyLong_AsLong(nbad) - lnsample = PyLong_AsLong(nsample) - if not PyErr_Occurred(): if lngood < 0: raise ValueError("ngood < 0") if lnbad < 0: @@ -5810,7 +5913,7 @@ cdef class _MKLRandomState: (( lnbad) != lnbad) or (( lnsample) != lnsample) ): - raise ValueError("All parameters should not exceed 2147483647") + raise ValueError("No parameters should exceed 2147483647") lntot = lngood + lnbad if lntot < lnsample: raise ValueError("ngood + nbad < nsample") @@ -5824,17 +5927,6 @@ cdef class _MKLRandomState: self.lock ) - PyErr_Clear() - - ongood = cnp.PyArray_FROM_OTF( - ngood, cnp.NPY_LONG, cnp.NPY_ARRAY_IN_ARRAY - ) - onbad = cnp.PyArray_FROM_OTF( - nbad, cnp.NPY_LONG, cnp.NPY_ARRAY_IN_ARRAY - ) - onsample = cnp.PyArray_FROM_OTF( - nsample, cnp.NPY_LONG, cnp.NPY_ARRAY_IN_ARRAY - ) if np.any(np.less(ongood, 0)): raise ValueError("ngood < 0") if np.any(np.less(onbad, 0)): @@ -5940,8 +6032,13 @@ cdef class _MKLRandomState: cdef cnp.ndarray op cdef double fp - fp = PyFloat_AsDouble(p) - if not PyErr_Occurred(): + op = cnp.PyArray_FROM_OTF( + p, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED + ) + + if cnp.PyArray_NDIM(op) == 0: + fp = PyFloat_AsDouble(p) + if fp <= 0.0: raise ValueError("p <= 0.0") if fp >= 1.0: @@ -5950,11 +6047,6 @@ cdef class _MKLRandomState: self.internal_state, irk_logseries_vec, size, fp, self.lock ) - PyErr_Clear() - - op = cnp.PyArray_FROM_OTF( - p, cnp.NPY_DOUBLE, cnp.NPY_ARRAY_ALIGNED - ) if np.any(np.less_equal(op, 0.0)): raise ValueError("p <= 0.0") if np.any(np.greater_equal(op, 1.0)): @@ -6400,17 +6492,20 @@ cdef class _MKLRandomState: cdef cnp.ndarray u "arrayObject_u" cdef double *u_data + if isinstance(x, np.ndarray) and not x.flags.writeable: + raise ValueError("array is read-only") + if (n == 0): return - u = self.random_sample(n-1) + u = self.random_sample(n - 1) u_data = cnp.PyArray_DATA(u) if type(x) is np.ndarray and x.ndim == 1 and x.size: # Fast, statically typed path: shuffle the underlying buffer. # Only for non-empty, 1d objects of class ndarray (subclasses such # as MaskedArrays may not support this approach). - x_ptr = x.ctypes.data + x_ptr = cnp.PyArray_BYTES(x) stride = x.strides[0] itemsize = x.dtype.itemsize # As the array x could contain python objects we use a buffer @@ -6418,7 +6513,7 @@ cdef class _MKLRandomState: # 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 = buf.ctypes.data + buf_ptr = cnp.PyArray_BYTES(buf) with self.lock: # We trick gcc into providing a specialized implementation for # the most common case, yielding a ~33% performance improvement. @@ -6431,9 +6526,19 @@ cdef class _MKLRandomState: self._shuffle_raw( n, itemsize, stride, x_ptr, buf_ptr, u_data ) - elif isinstance(x, np.ndarray) and x.ndim > 1 and x.size: - # Multidimensional ndarrays require a bounce buffer. - buf = np.empty_like(x[0]) + elif isinstance(x, np.ndarray): + if x.size == 0: + # shuffling is a no-op + return + + if x.ndim == 1 and x.dtype.type is np.object_: + warnings.warn( + "Shuffling a one dimensional array subclass containing " + "objects gives incorrect results for most array " + "subclasses.", + UserWarning, stacklevel=1) # Cython adds no stacklevel + + buf = np.empty_like(x[0, ...]) with self.lock: for i in reversed(range(1, n)): j = floor((i + 1) * u_data[i - 1]) @@ -6443,6 +6548,14 @@ cdef class _MKLRandomState: x[i] = buf else: # Untyped path. + if not isinstance(x, Sequence): + warnings.warn( + f"you are shuffling a '{type(x).__name__}' object " + "which is not a subclass of 'Sequence'; " + "`shuffle` is not guaranteed to behave correctly. " + "E.g., non-numpy array/tensor objects with view semantics " + "may contain duplicates after shuffling.", + UserWarning, stacklevel=1) # Cython does not add a level with self.lock: for i in reversed(range(1, n)): j = floor((i + 1) * u_data[i - 1]) @@ -6606,53 +6719,6 @@ cdef class MKLRandomState(_MKLRandomState): f"for {str(_brng_id_to_name(brng_id))}" ) - def tomaxint(self, size=None): - """ - tomaxint(size=None) - - Return a sample of uniformly distributed random integers in the interval - [0, ``np.iinfo("long").max``]. - - Parameters - ---------- - 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. Default is None, in which case a - single value is returned. - - Returns - ------- - out : ndarray - Drawn samples, with shape `size`. - - See Also - -------- - randint : Uniform sampling over a given half-open interval of integers. - random_integers : Uniform sampling over a given closed interval of - integers. - - Examples - -------- - >>> RS = mkl_random.MKLRandomState() # need a MKLRandomState object - >>> RS.tomaxint((2,2,2)) - array([[[1170048599, 1600360186], - [ 739731006, 1947757578]], - [[1871712945, 752307660], - [1601631370, 1479324245]]]) - >>> import sys - >>> sys.maxint - 2147483647 - >>> RS.tomaxint((2,2,2)) < sys.maxint - array([[[ True, True], - [ True, True]], - [[ True, True], - [ True, True]]], dtype=bool) - - """ - return vec_long_disc0_array( - self.internal_state, irk_long_vec, size, self.lock - ) - def randint_untyped(self, low, high=None, size=None): """ randint_untyped(low, high=None, size=None, dtype=int) diff --git a/mkl_random/tests/third_party/__init__.py b/mkl_random/tests/third_party/__init__.py new file mode 100644 index 0000000..ff61104 --- /dev/null +++ b/mkl_random/tests/third_party/__init__.py @@ -0,0 +1,24 @@ +# Copyright (c) 2017, Intel Corporation +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# * Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# * Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# * Neither the name of Intel Corporation nor the names of its contributors +# may be used to endorse or promote products derived from this software +# without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE +# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/mkl_random/tests/third_party/test_numpy_random.py b/mkl_random/tests/third_party/test_numpy_random.py new file mode 100644 index 0000000..71ec88e --- /dev/null +++ b/mkl_random/tests/third_party/test_numpy_random.py @@ -0,0 +1,1130 @@ +# This file includes tests from numpy.random module: +# https://github.com/numpy/numpy/blob/main/numpy/random/tests/test_random.py + +import warnings + +import numpy as np +import pytest +from numpy.testing import ( + assert_, + assert_array_equal, + assert_equal, + assert_no_warnings, + assert_raises, +) + +import mkl_random.interfaces.numpy_random as mkl_random + + +class TestSeed: + # MKL can't guarantee that seed will produce the same results for different + # architectures or platforms, so we test that the seed is accepted + # and produces a result + def test_scalar(self): + s = mkl_random.RandomState(0) + assert isinstance(s.randint(1000), int) + + def test_array(self): + s = mkl_random.RandomState(range(10)) + assert isinstance(s.randint(1000), int) + s = mkl_random.RandomState(np.arange(10)) + assert isinstance(s.randint(1000), int) + s = mkl_random.RandomState([0]) + assert isinstance(s.randint(1000), int) + s = mkl_random.RandomState([4294967295]) + assert isinstance(s.randint(1000), int) + + def test_invalid_scalar(self): + # seed must be an unsigned 32 bit integer + assert_raises(TypeError, mkl_random.RandomState, -0.5) + assert_raises(ValueError, mkl_random.RandomState, -1) + + def test_invalid_array(self): + # seed must be an unsigned 32 bit integer + assert_raises(TypeError, mkl_random.RandomState, [-0.5]) + assert_raises(ValueError, mkl_random.RandomState, [-1]) + assert_raises(ValueError, mkl_random.RandomState, [4294967296]) + assert_raises(ValueError, mkl_random.RandomState, [1, 2, 4294967296]) + assert_raises(ValueError, mkl_random.RandomState, [1, -2, 4294967296]) + + def test_invalid_array_shape(self): + # gh-9832 + assert_raises( + ValueError, mkl_random.RandomState, np.array([], dtype=np.int64) + ) + assert_raises(ValueError, mkl_random.RandomState, [[1, 2, 3]]) + assert_raises( + ValueError, mkl_random.RandomState, [[1, 2, 3], [4, 5, 6]] + ) + + +class TestBinomial: + def test_n_zero(self): + # Tests the corner case of n == 0 for the binomial distribution. + # binomial(0, p) should be zero for any p in [0, 1]. + # This test addresses issue #3480. + zeros = np.zeros(2, dtype="int") + for p in [0, 0.5, 1]: + assert_(mkl_random.binomial(0, p) == 0) + assert_array_equal(mkl_random.binomial(zeros, p), zeros) + + def test_p_is_nan(self): + # Issue #4571. + assert_raises(ValueError, mkl_random.binomial, 1, np.nan) + + +class TestMultinomial: + def test_basic(self): + mkl_random.multinomial(100, [0.2, 0.8]) + + def test_zero_probability(self): + mkl_random.multinomial(100, [0.2, 0.8, 0.0, 0.0, 0.0]) + + def test_int_negative_interval(self): + assert_(-5 <= mkl_random.randint(-5, -1) < -1) + x = mkl_random.randint(-5, -1, 5) + assert_(np.all(-5 <= x)) + assert_(np.all(x < -1)) + + def test_size(self): + # gh-3173 + p = [0.5, 0.5] + assert_equal(mkl_random.multinomial(1, p, np.uint32(1)).shape, (1, 2)) + assert_equal(mkl_random.multinomial(1, p, np.uint32(1)).shape, (1, 2)) + assert_equal(mkl_random.multinomial(1, p, np.uint32(1)).shape, (1, 2)) + assert_equal(mkl_random.multinomial(1, p, [2, 2]).shape, (2, 2, 2)) + assert_equal(mkl_random.multinomial(1, p, (2, 2)).shape, (2, 2, 2)) + assert_equal( + mkl_random.multinomial(1, p, np.array((2, 2))).shape, (2, 2, 2) + ) + + assert_raises(TypeError, mkl_random.multinomial, 1, p, float(1)) + + def test_multidimensional_pvals(self): + assert_raises(ValueError, mkl_random.multinomial, 10, [[0, 1]]) + assert_raises(ValueError, mkl_random.multinomial, 10, [[0], [1]]) + assert_raises( + ValueError, mkl_random.multinomial, 10, [[[0], [1]], [[1], [0]]] + ) + assert_raises( + ValueError, mkl_random.multinomial, 10, np.array([[0, 1], [1, 0]]) + ) + + +class TestSetState: + def _create_rng(self): + seed = 1234567890 + prng = mkl_random.RandomState(seed) + state = prng.get_state() + return prng, state + + def test_basic(self): + prng, state = self._create_rng() + old = prng.tomaxint(16) + prng.set_state(state) + new = prng.tomaxint(16) + assert_(np.all(old == new)) + + def test_negative_binomial(self): + # Ensure that the negative binomial results take floating point + # arguments without truncation. + prng, _ = self._create_rng() + prng.negative_binomial(0.5, 0.5) + + def test_set_invalid_state(self): + # gh-25402 + prng, _ = self._create_rng() + with pytest.raises(IndexError): + prng.set_state(()) + + +class TestRandint: + # valid integer/boolean types + itype = [ + np.bool, + np.int8, + np.uint8, + np.int16, + np.uint16, + np.int32, + np.uint32, + np.int64, + np.uint64, + ] + + def test_unsupported_type(self): + rng = mkl_random.RandomState() + assert_raises(TypeError, rng.randint, 1, dtype=float) + + def test_bounds_checking(self): + rng = mkl_random.RandomState() + for dt in self.itype: + 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 + assert_raises(ValueError, rng.randint, lbnd - 1, ubnd, dtype=dt) + assert_raises(ValueError, rng.randint, lbnd, ubnd + 1, dtype=dt) + assert_raises(ValueError, rng.randint, ubnd, lbnd, dtype=dt) + assert_raises(ValueError, rng.randint, 1, 0, dtype=dt) + + def test_rng_zero_and_extremes(self): + rng = mkl_random.RandomState() + for dt in self.itype: + 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(rng.randint(tgt, tgt + 1, size=1000, dtype=dt), tgt) + + tgt = lbnd + assert_equal(rng.randint(tgt, tgt + 1, size=1000, dtype=dt), tgt) + + tgt = (lbnd + ubnd) // 2 + assert_equal(rng.randint(tgt, tgt + 1, size=1000, dtype=dt), tgt) + + def test_full_range(self): + # Test for ticket #1690 + rng = mkl_random.RandomState() + + for dt in self.itype: + 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 + + try: + rng.randint(lbnd, ubnd, dtype=dt) + except Exception as e: + raise AssertionError( + "No error should have been raised, " + "but one was with the following " + f"message:\n\n{str(e)}" + ) + + def test_in_bounds_fuzz(self): + rng = mkl_random.RandomState() + + for dt in self.itype[1:]: + for ubnd in [4, 8, 16]: + vals = rng.randint(2, ubnd, size=2**16, dtype=dt) + assert_(vals.max() < ubnd) + assert_(vals.min() >= 2) + + vals = rng.randint(0, 2, size=2**16, dtype=np.bool) + + assert_(vals.max() < 2) + assert_(vals.min() >= 0) + + def test_int64_uint64_corner_case(self): + # When stored in Numpy arrays, `lbnd` is casted + # as np.int64, and `ubnd` is casted as np.uint64. + # Checking whether `lbnd` >= `ubnd` used to be + # done solely via direct comparison, which is incorrect + # because when Numpy tries to compare both numbers, + # it casts both to np.float64 because there is + # no integer superset of np.int64 and np.uint64. However, + # `ubnd` is too large to be represented in np.float64, + # causing it be round down to np.iinfo(np.int64).max, + # leading to a ValueError because `lbnd` now equals + # the new `ubnd`. + + dt = np.int64 + tgt = np.iinfo(np.int64).max + lbnd = np.int64(np.iinfo(np.int64).max) + ubnd = np.uint64(np.iinfo(np.int64).max + 1) + + # None of these function calls should + # generate a ValueError now. + actual = mkl_random.randint(lbnd, ubnd, dtype=dt) + assert_equal(actual, tgt) + + def test_respect_dtype_singleton(self): + # See gh-7203 + rng = mkl_random.RandomState() + for dt in self.itype: + 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 + + sample = rng.randint(lbnd, ubnd, dtype=dt) + assert_equal(sample.dtype, np.dtype(dt)) + + for dt in (bool, int): + # 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 = rng.randint(lbnd, ubnd, dtype=dt) + assert_(not hasattr(sample, "dtype")) + assert_equal(type(sample), dt) + + +class TestRandomDist: + 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 attempting + # to generate this integer. + with pytest.warns(DeprecationWarning): + actual = mkl_random.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(): + warnings.simplefilter("error", DeprecationWarning) + + # DeprecationWarning raised with high == None + assert_raises( + DeprecationWarning, + mkl_random.random_integers, + np.iinfo("l").max, + ) + + # DeprecationWarning raised with high != None + assert_raises( + DeprecationWarning, + mkl_random.random_integers, + np.iinfo("l").max, + np.iinfo("l").max, + ) + + def test_choice_exceptions(self): + sample = mkl_random.choice + assert_raises(ValueError, sample, -1, 3) + assert_raises(ValueError, sample, 3.0, 3) + assert_raises(ValueError, sample, [[1, 2], [3, 4]], 3) + assert_raises(ValueError, sample, [], 3) + assert_raises( + ValueError, sample, [1, 2, 3, 4], 3, p=[[0.25, 0.25], [0.25, 0.25]] + ) + assert_raises(ValueError, sample, [1, 2], 3, p=[0.4, 0.4, 0.2]) + assert_raises(ValueError, sample, [1, 2], 3, p=[1.1, -0.1]) + assert_raises(ValueError, sample, [1, 2], 3, p=[0.4, 0.4]) + assert_raises(ValueError, sample, [1, 2, 3], 4, replace=False) + # gh-13087 + assert_raises(ValueError, sample, [1, 2, 3], -2, replace=False) + assert_raises(ValueError, sample, [1, 2, 3], (-1,), replace=False) + assert_raises(ValueError, sample, [1, 2, 3], (-1, 1), replace=False) + assert_raises( + ValueError, sample, [1, 2, 3], 2, replace=False, p=[1, 0, 0] + ) + + def test_choice_return_shape(self): + p = [0.1, 0.9] + # Check scalar + assert_(np.isscalar(mkl_random.choice(2, replace=True))) + assert_(np.isscalar(mkl_random.choice(2, replace=False))) + assert_(np.isscalar(mkl_random.choice(2, replace=True, p=p))) + assert_(np.isscalar(mkl_random.choice(2, replace=False, p=p))) + assert_(np.isscalar(mkl_random.choice([1, 2], replace=True))) + assert_(mkl_random.choice([None], replace=True) is None) + a = np.array([1, 2]) + arr = np.empty(1, dtype=object) + arr[0] = a + assert_(mkl_random.choice(arr, replace=True) is a) + + # Check 0-d array + s = () + assert_(not np.isscalar(mkl_random.choice(2, s, replace=True))) + assert_(not np.isscalar(mkl_random.choice(2, s, replace=False))) + assert_(not np.isscalar(mkl_random.choice(2, s, replace=True, p=p))) + assert_(not np.isscalar(mkl_random.choice(2, s, replace=False, p=p))) + assert_(not np.isscalar(mkl_random.choice([1, 2], s, replace=True))) + assert_(mkl_random.choice([None], s, replace=True).ndim == 0) + a = np.array([1, 2]) + arr = np.empty(1, dtype=object) + arr[0] = a + assert_(mkl_random.choice(arr, s, replace=True).item() is a) + + # Check multi dimensional array + s = (2, 3) + p = [0.1, 0.1, 0.1, 0.1, 0.4, 0.2] + assert_equal(mkl_random.choice(6, s, replace=True).shape, s) + assert_equal(mkl_random.choice(6, s, replace=False).shape, s) + assert_equal(mkl_random.choice(6, s, replace=True, p=p).shape, s) + assert_equal(mkl_random.choice(6, s, replace=False, p=p).shape, s) + assert_equal(mkl_random.choice(np.arange(6), s, replace=True).shape, s) + + # Check zero-size + assert_equal(mkl_random.randint(0, 0, size=(3, 0, 4)).shape, (3, 0, 4)) + assert_equal(mkl_random.randint(0, -10, size=0).shape, (0,)) + assert_equal(mkl_random.randint(10, 10, size=0).shape, (0,)) + assert_equal(mkl_random.choice(0, size=0).shape, (0,)) + assert_equal(mkl_random.choice([], size=(0,)).shape, (0,)) + assert_equal( + mkl_random.choice(["a", "b"], size=(3, 0, 4)).shape, (3, 0, 4) + ) + assert_raises(ValueError, mkl_random.choice, [], 10) + + def test_choice_nan_probabilities(self): + a = np.array([42, 1, 2]) + p = [None, None, None] + assert_raises(ValueError, mkl_random.choice, a, p=p) + + def test_shuffle(self): + # Test lists, arrays (of various dtypes), and multidimensional versions + # of both, c-contiguous or not: + for conv in [ + lambda x: np.array([]), + lambda x: x, + lambda x: np.asarray(x).astype(np.int8), + lambda x: np.asarray(x).astype(np.float32), + lambda x: np.asarray(x).astype(np.complex64), + lambda x: np.asarray(x).astype(object), + lambda x: [(i, i) for i in x], + lambda x: np.asarray([[i, i] for i in x]), + lambda x: np.vstack([x, x]).T, + # gh-11442 + lambda x: ( + np.asarray([(i, i) for i in x], [("a", int), ("b", int)]).view( + np.recarray + ) + ), + # gh-4270 + lambda x: np.asarray( + [(i, i) for i in x], [("a", object), ("b", np.int32)] + ), + ]: + rng = mkl_random.RandomState() + alist = conv([1, 2, 3, 4, 5, 6, 7, 8, 9, 0]) + # Do not validate against expected results as we cannot guarantee + # consistency across platforms or architectures. + # This test is just to check that it runs on all types + rng.shuffle(alist) + + @pytest.mark.parametrize("random", [mkl_random, mkl_random.RandomState()]) + def test_shuffle_untyped_warning(self, random): + # Create a dict works like a sequence but isn't one + values = {0: 0, 1: 1, 2: 2, 3: 3, 4: 4, 5: 5, 6: 6} + with pytest.warns( + UserWarning, match="you are shuffling a 'dict' object" + ): + random.shuffle(values) + + @pytest.mark.parametrize("random", [mkl_random, mkl_random.RandomState()]) + @pytest.mark.parametrize("use_array_like", [True, False]) + def test_shuffle_no_object_unpacking(self, random, use_array_like): + class MyArr(np.ndarray): + pass + + items = [ + None, + np.array([3]), + np.float64(3), + np.array(10), + np.float64(7), + ] + arr = np.array(items, dtype=object) + item_ids = {id(i) for i in items} + if use_array_like: + arr = arr.view(MyArr) + + # The array was created fine, and did not modify any objects: + assert all(id(i) in item_ids for i in arr) + + if use_array_like: + with pytest.warns( + UserWarning, match="Shuffling a one dimensional array.*" + ): + random.shuffle(arr) + else: + random.shuffle(arr) + assert all(id(i) in item_ids for i in arr) + + def test_shuffle_memoryview(self): + # gh-18273 + # allow graceful handling of memoryviews + # (treat the same as arrays) + rng = mkl_random.RandomState() + a = np.arange(5).data + rng.shuffle(a) + + def test_shuffle_not_writeable(self): + a = np.zeros(3) + a.flags.writeable = False + with pytest.raises(ValueError, match="read-only"): + mkl_random.shuffle(a) + + def test_dirichlet_size(self): + # gh-3173 + p = np.array([51.72840233779265162, 39.74494232180943953]) + assert_equal(mkl_random.dirichlet(p, np.uint32(1)).shape, (1, 2)) + assert_equal(mkl_random.dirichlet(p, np.uint32(1)).shape, (1, 2)) + assert_equal(mkl_random.dirichlet(p, np.uint32(1)).shape, (1, 2)) + assert_equal(mkl_random.dirichlet(p, [2, 2]).shape, (2, 2, 2)) + assert_equal(mkl_random.dirichlet(p, (2, 2)).shape, (2, 2, 2)) + assert_equal(mkl_random.dirichlet(p, np.array((2, 2))).shape, (2, 2, 2)) + + assert_raises(TypeError, mkl_random.dirichlet, p, float(1)) + + def test_dirichlet_bad_alpha(self): + # gh-2089 + alpha = np.array([5.4e-01, -1.0e-16]) + assert_raises(ValueError, mkl_random.dirichlet, alpha) + + # gh-15876 + assert_raises(ValueError, mkl_random.dirichlet, [[5, 1]]) + assert_raises(ValueError, mkl_random.dirichlet, [[5], [1]]) + assert_raises( + ValueError, mkl_random.dirichlet, [[[5], [1]], [[1], [5]]] + ) + assert_raises( + ValueError, mkl_random.dirichlet, np.array([[5, 1], [1, 5]]) + ) + + def test_multivariate_normal_warnings(self): + rng = mkl_random.RandomState() + + # Check that non positive-semidefinite covariance warns with + # RuntimeWarning + mean = [0, 0] + cov = [[1, 2], [2, 1]] + pytest.warns(RuntimeWarning, rng.multivariate_normal, mean, cov) + + # and that it doesn't warn with RuntimeWarning check_valid='ignore' + assert_no_warnings( + rng.multivariate_normal, mean, cov, check_valid="ignore" + ) + + # and that it raises with RuntimeWarning check_valid='raises' + assert_raises( + ValueError, rng.multivariate_normal, mean, cov, check_valid="raise" + ) + + cov = np.array([[1, 0.1], [0.1, 1]], dtype=np.float32) + with warnings.catch_warnings(): + warnings.simplefilter("error") + rng.multivariate_normal(mean, cov) + + def test_poisson_exceptions(self): + lambig = np.iinfo("l").max + lamneg = -1 + assert_raises(ValueError, mkl_random.poisson, lamneg) + assert_raises(ValueError, mkl_random.poisson, [lamneg] * 10) + assert_raises(ValueError, mkl_random.poisson, lambig) + assert_raises(ValueError, mkl_random.poisson, [lambig] * 10) + + # TODO: revisit test after experimenting with range calculation in uniform + @pytest.mark.skip("Uniform does not overflow identically to NumPy") + def test_uniform_range_bounds(self): + fmin = np.finfo("float").min + fmax = np.finfo("float").max + + func = mkl_random.uniform + assert_raises(OverflowError, func, -np.inf, 0) + assert_raises(OverflowError, func, 0, np.inf) + assert_raises(OverflowError, func, fmin, fmax) + assert_raises(OverflowError, func, [-np.inf], [0]) + assert_raises(OverflowError, func, [0], [np.inf]) + + # (fmax / 1e17) - fmin is within range, so this should not throw + # account for i386 extended precision DBL_MAX / 1e17 + DBL_MAX > + # DBL_MAX by increasing fmin a bit + mkl_random.uniform(low=np.nextafter(fmin, 1), high=fmax / 1e17) + + # TODO: revisit after changing conversion logic + @pytest.mark.skip("mkl_random casts via NumPy instead of throwing") + def test_scalar_exception_propagation(self): + # Tests that exceptions are correctly propagated in distributions + # when called with objects that throw exceptions when converted to + # scalars. + # + # Regression test for gh: 8865 + + class ThrowingFloat(np.ndarray): + def __float__(self): + raise TypeError + + throwing_float = np.array(1.0).view(ThrowingFloat) + assert_raises( + TypeError, mkl_random.uniform, throwing_float, throwing_float + ) + + class ThrowingInteger(np.ndarray): + def __int__(self): + raise TypeError + + __index__ = __int__ + + throwing_int = np.array(1).view(ThrowingInteger) + assert_raises(TypeError, mkl_random.hypergeometric, throwing_int, 1, 1) + + def test_vonmises_small(self): + # check infinite loop, gh-4720 + mkl_random.seed() + r = mkl_random.vonmises(mu=0.0, kappa=1.1e-8, size=10**6) + np.testing.assert_(np.isfinite(r).all()) + + +class TestBroadcast: + def test_uniform(self): + low = [0] + high = [1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.uniform(low * 3, high).shape, desired_shape) + assert_equal(rng.uniform(low, high * 3).shape, desired_shape) + + def test_normal(self): + loc = [0] + scale = [1] + bad_scale = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.normal(loc * 3, scale).shape, desired_shape) + assert_raises(ValueError, rng.normal, loc * 3, bad_scale) + + assert_equal(rng.normal(loc, scale * 3).shape, desired_shape) + assert_raises(ValueError, rng.normal, loc, bad_scale * 3) + + def test_beta(self): + a = [1] + b = [2] + bad_a = [-1] + bad_b = [-2] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.beta(a * 3, b).shape, desired_shape) + assert_raises(ValueError, rng.beta, bad_a * 3, b) + assert_raises(ValueError, rng.beta, a * 3, bad_b) + + assert_equal(rng.beta(a, b * 3).shape, desired_shape) + assert_raises(ValueError, rng.beta, bad_a, b * 3) + assert_raises(ValueError, rng.beta, a, bad_b * 3) + + def test_exponential(self): + scale = [1] + bad_scale = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.exponential(scale * 3).shape, desired_shape) + assert_raises(ValueError, rng.exponential, bad_scale * 3) + + def test_standard_gamma(self): + shape = [1] + bad_shape = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.standard_gamma(shape * 3).shape, desired_shape) + assert_raises(ValueError, rng.standard_gamma, bad_shape * 3) + + def test_gamma(self): + shape = [1] + scale = [2] + bad_shape = [-1] + bad_scale = [-2] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.gamma(shape * 3, scale).shape, desired_shape) + assert_raises(ValueError, rng.gamma, bad_shape * 3, scale) + assert_raises(ValueError, rng.gamma, shape * 3, bad_scale) + + assert_equal(rng.gamma(shape, scale * 3).shape, desired_shape) + assert_raises(ValueError, rng.gamma, bad_shape, scale * 3) + assert_raises(ValueError, rng.gamma, shape, bad_scale * 3) + + def test_f(self): + dfnum = [1] + dfden = [2] + bad_dfnum = [-1] + bad_dfden = [-2] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.f(dfnum * 3, dfden).shape, desired_shape) + assert_raises(ValueError, rng.f, bad_dfnum * 3, dfden) + assert_raises(ValueError, rng.f, dfnum * 3, bad_dfden) + + assert_equal(rng.f(dfnum, dfden * 3).shape, desired_shape) + assert_raises(ValueError, rng.f, bad_dfnum, dfden * 3) + assert_raises(ValueError, rng.f, dfnum, bad_dfden * 3) + + def test_noncentral_f(self): + dfnum = [2] + dfden = [3] + nonc = [4] + bad_dfnum = [0] + bad_dfden = [-1] + bad_nonc = [-2] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal( + rng.noncentral_f(dfnum * 3, dfden, nonc).shape, desired_shape + ) + assert_raises(ValueError, rng.noncentral_f, bad_dfnum * 3, dfden, nonc) + assert_raises(ValueError, rng.noncentral_f, dfnum * 3, bad_dfden, nonc) + assert_raises(ValueError, rng.noncentral_f, dfnum * 3, dfden, bad_nonc) + + assert_equal( + rng.noncentral_f(dfnum, dfden * 3, nonc).shape, desired_shape + ) + assert_raises(ValueError, rng.noncentral_f, bad_dfnum, dfden * 3, nonc) + assert_raises(ValueError, rng.noncentral_f, dfnum, bad_dfden * 3, nonc) + assert_raises(ValueError, rng.noncentral_f, dfnum, dfden * 3, bad_nonc) + + assert_equal( + rng.noncentral_f(dfnum, dfden, nonc * 3).shape, desired_shape + ) + assert_raises(ValueError, rng.noncentral_f, bad_dfnum, dfden, nonc * 3) + assert_raises(ValueError, rng.noncentral_f, dfnum, bad_dfden, nonc * 3) + assert_raises(ValueError, rng.noncentral_f, dfnum, dfden, bad_nonc * 3) + + def test_chisquare(self): + df = [1] + bad_df = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.chisquare(df * 3).shape, desired_shape) + assert_raises(ValueError, rng.chisquare, bad_df * 3) + + def test_noncentral_chisquare(self): + df = [1] + nonc = [2] + bad_df = [-1] + bad_nonc = [-2] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal( + rng.noncentral_chisquare(df * 3, nonc).shape, desired_shape + ) + assert_raises(ValueError, rng.noncentral_chisquare, bad_df * 3, nonc) + assert_raises(ValueError, rng.noncentral_chisquare, df * 3, bad_nonc) + + assert_equal( + rng.noncentral_chisquare(df, nonc * 3).shape, desired_shape + ) + assert_raises(ValueError, rng.noncentral_chisquare, bad_df, nonc * 3) + assert_raises(ValueError, rng.noncentral_chisquare, df, bad_nonc * 3) + + def test_standard_t(self): + df = [1] + bad_df = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.standard_t(df * 3).shape, desired_shape) + assert_raises(ValueError, rng.standard_t, bad_df * 3) + + def test_vonmises(self): + mu = [2] + kappa = [1] + bad_kappa = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.vonmises(mu * 3, kappa).shape, desired_shape) + assert_raises(ValueError, rng.vonmises, mu * 3, bad_kappa) + + assert_equal(rng.vonmises(mu, kappa * 3).shape, desired_shape) + assert_raises(ValueError, rng.vonmises, mu, bad_kappa * 3) + + def test_pareto(self): + a = [1] + bad_a = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.pareto(a * 3).shape, desired_shape) + assert_raises(ValueError, rng.pareto, bad_a * 3) + + def test_weibull(self): + a = [1] + bad_a = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.weibull(a * 3).shape, desired_shape) + assert_raises(ValueError, rng.weibull, bad_a * 3) + + def test_power(self): + a = [1] + bad_a = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.power(a * 3).shape, desired_shape) + assert_raises(ValueError, rng.power, bad_a * 3) + + def test_laplace(self): + loc = [0] + scale = [1] + bad_scale = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.laplace(loc * 3, scale).shape, desired_shape) + assert_raises(ValueError, rng.laplace, loc * 3, bad_scale) + + assert_equal(rng.laplace(loc, scale * 3).shape, desired_shape) + assert_raises(ValueError, rng.laplace, loc, bad_scale * 3) + + def test_gumbel(self): + loc = [0] + scale = [1] + bad_scale = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.gumbel(loc * 3, scale).shape, desired_shape) + assert_raises(ValueError, rng.gumbel, loc * 3, bad_scale) + + assert_equal(rng.gumbel(loc, scale * 3).shape, desired_shape) + assert_raises(ValueError, rng.gumbel, loc, bad_scale * 3) + + def test_logistic(self): + loc = [0] + scale = [1] + bad_scale = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.logistic(loc * 3, scale).shape, desired_shape) + assert_raises(ValueError, rng.logistic, loc * 3, bad_scale) + + assert_equal(rng.logistic(loc, scale * 3).shape, desired_shape) + assert_raises(ValueError, rng.logistic, loc, bad_scale * 3) + + def test_lognormal(self): + mean = [0] + sigma = [1] + bad_sigma = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.lognormal(mean * 3, sigma).shape, desired_shape) + assert_raises(ValueError, rng.lognormal, mean * 3, bad_sigma) + + assert_equal(rng.lognormal(mean, sigma * 3).shape, desired_shape) + assert_raises(ValueError, rng.lognormal, mean, bad_sigma * 3) + + def test_rayleigh(self): + scale = [1] + bad_scale = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.rayleigh(scale * 3).shape, desired_shape) + assert_raises(ValueError, rng.rayleigh, bad_scale * 3) + + def test_wald(self): + mean = [0.5] + scale = [1] + bad_mean = [0] + bad_scale = [-2] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.wald(mean * 3, scale).shape, desired_shape) + assert_raises(ValueError, rng.wald, bad_mean * 3, scale) + assert_raises(ValueError, rng.wald, mean * 3, bad_scale) + + assert_equal(rng.wald(mean, scale * 3).shape, desired_shape) + assert_raises(ValueError, rng.wald, bad_mean, scale * 3) + assert_raises(ValueError, rng.wald, mean, bad_scale * 3) + assert_raises(ValueError, rng.wald, 0.0, 1) + assert_raises(ValueError, rng.wald, 0.5, 0.0) + + def test_triangular(self): + left = [1] + right = [3] + mode = [2] + bad_left_one = [3] + bad_mode_one = [4] + bad_left_two, bad_mode_two = right * 2 + + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.triangular(left * 3, mode, right).shape, desired_shape) + assert_raises(ValueError, rng.triangular, bad_left_one * 3, mode, right) + assert_raises(ValueError, rng.triangular, left * 3, bad_mode_one, right) + assert_raises( + ValueError, rng.triangular, bad_left_two * 3, bad_mode_two, right + ) + + assert_equal(rng.triangular(left, mode * 3, right).shape, desired_shape) + assert_raises(ValueError, rng.triangular, bad_left_one, mode * 3, right) + assert_raises(ValueError, rng.triangular, left, bad_mode_one * 3, right) + assert_raises( + ValueError, rng.triangular, bad_left_two, bad_mode_two * 3, right + ) + + assert_equal(rng.triangular(left, mode, right * 3).shape, desired_shape) + assert_raises(ValueError, rng.triangular, bad_left_one, mode, right * 3) + assert_raises(ValueError, rng.triangular, left, bad_mode_one, right * 3) + assert_raises( + ValueError, rng.triangular, bad_left_two, bad_mode_two, right * 3 + ) + + def test_binomial(self): + n = [1] + p = [0.5] + bad_n = [-1] + bad_p_one = [-1] + bad_p_two = [1.5] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.binomial(n * 3, p).shape, desired_shape) + assert_raises(ValueError, rng.binomial, bad_n * 3, p) + assert_raises(ValueError, rng.binomial, n * 3, bad_p_one) + assert_raises(ValueError, rng.binomial, n * 3, bad_p_two) + + assert_equal(rng.binomial(n, p * 3).shape, desired_shape) + assert_raises(ValueError, rng.binomial, bad_n, p * 3) + assert_raises(ValueError, rng.binomial, n, bad_p_one * 3) + assert_raises(ValueError, rng.binomial, n, bad_p_two * 3) + + def test_negative_binomial(self): + n = [1] + p = [0.5] + bad_n = [-1] + bad_p_one = [-1] + bad_p_two = [1.5] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.negative_binomial(n * 3, p).shape, desired_shape) + assert_raises(ValueError, rng.negative_binomial, bad_n * 3, p) + assert_raises(ValueError, rng.negative_binomial, n * 3, bad_p_one) + assert_raises(ValueError, rng.negative_binomial, n * 3, bad_p_two) + + assert_equal(rng.negative_binomial(n, p * 3).shape, desired_shape) + assert_raises(ValueError, rng.negative_binomial, bad_n, p * 3) + assert_raises(ValueError, rng.negative_binomial, n, bad_p_one * 3) + assert_raises(ValueError, rng.negative_binomial, n, bad_p_two * 3) + + def test_poisson(self): + max_lam = mkl_random.RandomState()._poisson_lam_max + + lam = [1] + bad_lam_one = [-1] + bad_lam_two = [max_lam * 2] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.poisson(lam * 3).shape, desired_shape) + assert_raises(ValueError, rng.poisson, bad_lam_one * 3) + assert_raises(ValueError, rng.poisson, bad_lam_two * 3) + + def test_zipf(self): + a = [2] + bad_a = [0] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.zipf(a * 3).shape, desired_shape) + assert_raises(ValueError, rng.zipf, bad_a * 3) + with np.errstate(invalid="ignore"): + assert_raises(ValueError, rng.zipf, np.nan) + assert_raises(ValueError, rng.zipf, [0, 0, np.nan]) + + def test_geometric(self): + p = [0.5] + bad_p_one = [-1] + bad_p_two = [1.5] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.geometric(p * 3).shape, desired_shape) + assert_raises(ValueError, rng.geometric, bad_p_one * 3) + assert_raises(ValueError, rng.geometric, bad_p_two * 3) + + def test_hypergeometric(self): + ngood = [1] + nbad = [2] + nsample = [2] + bad_ngood = [-1] + bad_nbad = [-2] + bad_nsample_one = [0] + bad_nsample_two = [4] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal( + rng.hypergeometric(ngood * 3, nbad, nsample).shape, desired_shape + ) + assert_raises( + ValueError, rng.hypergeometric, bad_ngood * 3, nbad, nsample + ) + assert_raises( + ValueError, rng.hypergeometric, ngood * 3, bad_nbad, nsample + ) + assert_raises( + ValueError, rng.hypergeometric, ngood * 3, nbad, bad_nsample_one + ) + assert_raises( + ValueError, rng.hypergeometric, ngood * 3, nbad, bad_nsample_two + ) + + assert_equal( + rng.hypergeometric(ngood, nbad * 3, nsample).shape, desired_shape + ) + assert_raises( + ValueError, rng.hypergeometric, bad_ngood, nbad * 3, nsample + ) + assert_raises( + ValueError, rng.hypergeometric, ngood, bad_nbad * 3, nsample + ) + assert_raises( + ValueError, rng.hypergeometric, ngood, nbad * 3, bad_nsample_one + ) + assert_raises( + ValueError, rng.hypergeometric, ngood, nbad * 3, bad_nsample_two + ) + + assert_equal( + rng.hypergeometric(ngood, nbad, nsample * 3).shape, desired_shape + ) + assert_raises( + ValueError, rng.hypergeometric, bad_ngood, nbad, nsample * 3 + ) + assert_raises( + ValueError, rng.hypergeometric, ngood, bad_nbad, nsample * 3 + ) + assert_raises( + ValueError, rng.hypergeometric, ngood, nbad, bad_nsample_one * 3 + ) + assert_raises( + ValueError, rng.hypergeometric, ngood, nbad, bad_nsample_two * 3 + ) + + def test_logseries(self): + p = [0.5] + bad_p_one = [2] + bad_p_two = [-1] + desired_shape = (3,) + + rng = mkl_random.RandomState() + assert_equal(rng.logseries(p * 3).shape, desired_shape) + assert_raises(ValueError, rng.logseries, bad_p_one * 3) + assert_raises(ValueError, rng.logseries, bad_p_two * 3) + + +# See Issue #4263 +class TestSingleEltArrayInput: + def _create_arrays(self): + return np.array([2]), np.array([3]), np.array([4]), (1,) + + def test_one_arg_funcs(self): + argOne, _, _, tgtShape = self._create_arrays() + funcs = ( + mkl_random.exponential, + mkl_random.standard_gamma, + mkl_random.chisquare, + mkl_random.standard_t, + mkl_random.pareto, + mkl_random.weibull, + mkl_random.power, + mkl_random.rayleigh, + mkl_random.poisson, + mkl_random.zipf, + mkl_random.geometric, + mkl_random.logseries, + ) + + probfuncs = (mkl_random.geometric, mkl_random.logseries) + + for func in funcs: + if func in probfuncs: # p < 1.0 + out = func(np.array([0.5])) + + else: + out = func(argOne) + + assert_equal(out.shape, tgtShape) + + def test_two_arg_funcs(self): + argOne, argTwo, _, tgtShape = self._create_arrays() + funcs = ( + mkl_random.uniform, + mkl_random.normal, + mkl_random.beta, + mkl_random.gamma, + mkl_random.f, + mkl_random.noncentral_chisquare, + mkl_random.vonmises, + mkl_random.laplace, + mkl_random.gumbel, + mkl_random.logistic, + mkl_random.lognormal, + mkl_random.wald, + mkl_random.binomial, + mkl_random.negative_binomial, + ) + + probfuncs = (mkl_random.binomial, mkl_random.negative_binomial) + + for func in funcs: + if func in probfuncs: # p <= 1 + argTwo = np.array([0.5]) + + else: + argTwo = argTwo + + out = func(argOne, argTwo) + assert_equal(out.shape, tgtShape) + + out = func(argOne[0], argTwo) + assert_equal(out.shape, tgtShape) + + out = func(argOne, argTwo[0]) + assert_equal(out.shape, tgtShape) + + # TODO: fix randint to handle single arrays correctly, remove skip + @pytest.mark.skip("randint does not work with arrays") + def test_randint(self): + _, _, _, tgtShape = self._create_arrays() + itype = [ + bool, + np.int8, + np.uint8, + np.int16, + np.uint16, + np.int32, + np.uint32, + np.int64, + np.uint64, + ] + func = mkl_random.randint + high = np.array([1]) + low = np.array([0]) + + for dt in itype: + out = func(low, high, dtype=dt) + assert_equal(out.shape, tgtShape) + + out = func(low[0], high, dtype=dt) + assert_equal(out.shape, tgtShape) + + out = func(low, high[0], dtype=dt) + assert_equal(out.shape, tgtShape) + + def test_three_arg_funcs(self): + argOne, argTwo, argThree, tgtShape = self._create_arrays() + funcs = [ + mkl_random.noncentral_f, + mkl_random.triangular, + mkl_random.hypergeometric, + ] + + for func in funcs: + out = func(argOne, argTwo, argThree) + assert_equal(out.shape, tgtShape) + + out = func(argOne[0], argTwo, argThree) + assert_equal(out.shape, tgtShape) + + out = func(argOne, argTwo[0], argThree) + assert_equal(out.shape, tgtShape) diff --git a/pyproject.toml b/pyproject.toml index 6c136fc..145723c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -101,12 +101,10 @@ generated-members = ["RandomState", "min", "max"] [tool.setuptools] include-package-data = true +packages = ["mkl_random", "mkl_random.interfaces"] [tool.setuptools.dynamic] version = {attr = "mkl_random._version.__version__"} [tool.setuptools.package-data] -"mkl_random" = ["tests/*.py"] - -[tool.setuptools.packages.find] -include = ["mkl_random", "mkl_random.interfaces"] +"mkl_random" = ["tests/**/*.py"]