diff options
Diffstat (limited to 'numpy/random/mtrand/mtrand.pyx')
| -rw-r--r-- | numpy/random/mtrand/mtrand.pyx | 933 |
1 files changed, 542 insertions, 391 deletions
diff --git a/numpy/random/mtrand/mtrand.pyx b/numpy/random/mtrand/mtrand.pyx index 84b52079d..55138cba7 100644 --- a/numpy/random/mtrand/mtrand.pyx +++ b/numpy/random/mtrand/mtrand.pyx @@ -55,65 +55,66 @@ cdef extern from "randomkit.h": void rk_seed(unsigned long seed, rk_state *state) rk_error rk_randomseed(rk_state *state) unsigned long rk_random(rk_state *state) - long rk_long(rk_state *state) - unsigned long rk_ulong(rk_state *state) - unsigned long rk_interval(unsigned long max, rk_state *state) - double rk_double(rk_state *state) - void rk_fill(void *buffer, size_t size, rk_state *state) + long rk_long(rk_state *state) nogil + unsigned long rk_ulong(rk_state *state) nogil + unsigned long rk_interval(unsigned long max, rk_state *state) nogil + double rk_double(rk_state *state) nogil + void rk_fill(void *buffer, size_t size, rk_state *state) nogil rk_error rk_devfill(void *buffer, size_t size, int strong) rk_error rk_altfill(void *buffer, size_t size, int strong, - rk_state *state) - double rk_gauss(rk_state *state) + rk_state *state) nogil + double rk_gauss(rk_state *state) nogil cdef extern from "distributions.h": - - double rk_normal(rk_state *state, double loc, double scale) - double rk_standard_exponential(rk_state *state) - double rk_exponential(rk_state *state, double scale) - double rk_uniform(rk_state *state, double loc, double scale) - double rk_standard_gamma(rk_state *state, double shape) - double rk_gamma(rk_state *state, double shape, double scale) - double rk_beta(rk_state *state, double a, double b) - double rk_chisquare(rk_state *state, double df) - double rk_noncentral_chisquare(rk_state *state, double df, double nonc) - double rk_f(rk_state *state, double dfnum, double dfden) - double rk_noncentral_f(rk_state *state, double dfnum, double dfden, double nonc) - double rk_standard_cauchy(rk_state *state) - double rk_standard_t(rk_state *state, double df) - double rk_vonmises(rk_state *state, double mu, double kappa) - double rk_pareto(rk_state *state, double a) - double rk_weibull(rk_state *state, double a) - double rk_power(rk_state *state, double a) - double rk_laplace(rk_state *state, double loc, double scale) - double rk_gumbel(rk_state *state, double loc, double scale) - double rk_logistic(rk_state *state, double loc, double scale) - double rk_lognormal(rk_state *state, double mode, double sigma) - double rk_rayleigh(rk_state *state, double mode) - double rk_wald(rk_state *state, double mean, double scale) - double rk_triangular(rk_state *state, double left, double mode, double right) - - long rk_binomial(rk_state *state, long n, double p) - long rk_binomial_btpe(rk_state *state, long n, double p) - long rk_binomial_inversion(rk_state *state, long n, double p) - long rk_negative_binomial(rk_state *state, double n, double p) - long rk_poisson(rk_state *state, double lam) - long rk_poisson_mult(rk_state *state, double lam) - long rk_poisson_ptrs(rk_state *state, double lam) - long rk_zipf(rk_state *state, double a) - long rk_geometric(rk_state *state, double p) - long rk_hypergeometric(rk_state *state, long good, long bad, long sample) - long rk_logseries(rk_state *state, double p) - -ctypedef double (* rk_cont0)(rk_state *state) -ctypedef double (* rk_cont1)(rk_state *state, double a) -ctypedef double (* rk_cont2)(rk_state *state, double a, double b) -ctypedef double (* rk_cont3)(rk_state *state, double a, double b, double c) - -ctypedef long (* rk_disc0)(rk_state *state) -ctypedef long (* rk_discnp)(rk_state *state, long n, double p) -ctypedef long (* rk_discdd)(rk_state *state, double n, double p) -ctypedef long (* rk_discnmN)(rk_state *state, long n, long m, long N) -ctypedef long (* rk_discd)(rk_state *state, double a) + # do not need the GIL, but they do need a lock on the state !! */ + + double rk_normal(rk_state *state, double loc, double scale) nogil + double rk_standard_exponential(rk_state *state) nogil + double rk_exponential(rk_state *state, double scale) nogil + double rk_uniform(rk_state *state, double loc, double scale) nogil + double rk_standard_gamma(rk_state *state, double shape) nogil + double rk_gamma(rk_state *state, double shape, double scale) nogil + double rk_beta(rk_state *state, double a, double b) nogil + double rk_chisquare(rk_state *state, double df) nogil + double rk_noncentral_chisquare(rk_state *state, double df, double nonc) nogil + double rk_f(rk_state *state, double dfnum, double dfden) nogil + double rk_noncentral_f(rk_state *state, double dfnum, double dfden, double nonc) nogil + double rk_standard_cauchy(rk_state *state) nogil + double rk_standard_t(rk_state *state, double df) nogil + double rk_vonmises(rk_state *state, double mu, double kappa) nogil + double rk_pareto(rk_state *state, double a) nogil + double rk_weibull(rk_state *state, double a) nogil + double rk_power(rk_state *state, double a) nogil + double rk_laplace(rk_state *state, double loc, double scale) nogil + double rk_gumbel(rk_state *state, double loc, double scale) nogil + double rk_logistic(rk_state *state, double loc, double scale) nogil + double rk_lognormal(rk_state *state, double mode, double sigma) nogil + double rk_rayleigh(rk_state *state, double mode) nogil + double rk_wald(rk_state *state, double mean, double scale) nogil + double rk_triangular(rk_state *state, double left, double mode, double right) nogil + + long rk_binomial(rk_state *state, long n, double p) nogil + long rk_binomial_btpe(rk_state *state, long n, double p) nogil + long rk_binomial_inversion(rk_state *state, long n, double p) nogil + long rk_negative_binomial(rk_state *state, double n, double p) nogil + long rk_poisson(rk_state *state, double lam) nogil + long rk_poisson_mult(rk_state *state, double lam) nogil + long rk_poisson_ptrs(rk_state *state, double lam) nogil + long rk_zipf(rk_state *state, double a) nogil + long rk_geometric(rk_state *state, double p) nogil + long rk_hypergeometric(rk_state *state, long good, long bad, long sample) nogil + long rk_logseries(rk_state *state, double p) nogil + +ctypedef double (* rk_cont0)(rk_state *state) nogil +ctypedef double (* rk_cont1)(rk_state *state, double a) nogil +ctypedef double (* rk_cont2)(rk_state *state, double a, double b) nogil +ctypedef double (* rk_cont3)(rk_state *state, double a, double b, double c) nogil + +ctypedef long (* rk_disc0)(rk_state *state) nogil +ctypedef long (* rk_discnp)(rk_state *state, long n, double p) nogil +ctypedef long (* rk_discdd)(rk_state *state, double n, double p) nogil +ctypedef long (* rk_discnmN)(rk_state *state, long n, long m, long N) nogil +ctypedef long (* rk_discd)(rk_state *state, double a) nogil cdef extern from "initarray.h": @@ -125,8 +126,11 @@ import_array() import numpy as np import operator +import warnings +from threading import Lock -cdef object cont0_array(rk_state *state, rk_cont0 func, object size): +cdef object cont0_array(rk_state *state, rk_cont0 func, object size, + object lock): cdef double *array_data cdef ndarray array "arrayObject" cdef npy_intp length @@ -138,12 +142,14 @@ cdef object cont0_array(rk_state *state, rk_cont0 func, object size): array = <ndarray>np.empty(size, np.float64) length = PyArray_SIZE(array) array_data = <double *>PyArray_DATA(array) - for i from 0 <= i < length: - array_data[i] = func(state) + with lock, nogil: + for i from 0 <= i < length: + array_data[i] = func(state) return array -cdef object cont1_array_sc(rk_state *state, rk_cont1 func, object size, double a): +cdef object cont1_array_sc(rk_state *state, rk_cont1 func, object size, double a, + object lock): cdef double *array_data cdef ndarray array "arrayObject" cdef npy_intp length @@ -155,11 +161,13 @@ cdef object cont1_array_sc(rk_state *state, rk_cont1 func, object size, double a array = <ndarray>np.empty(size, np.float64) length = PyArray_SIZE(array) array_data = <double *>PyArray_DATA(array) - for i from 0 <= i < length: - array_data[i] = func(state, a) + with lock, nogil: + for i from 0 <= i < length: + array_data[i] = func(state, a) return array -cdef object cont1_array(rk_state *state, rk_cont1 func, object size, ndarray oa): +cdef object cont1_array(rk_state *state, rk_cont1 func, object size, + ndarray oa, object lock): cdef double *array_data cdef double *oa_data cdef ndarray array "arrayObject" @@ -174,9 +182,10 @@ cdef object cont1_array(rk_state *state, rk_cont1 func, object size, ndarray oa) length = PyArray_SIZE(array) array_data = <double *>PyArray_DATA(array) itera = <flatiter>PyArray_IterNew(<object>oa) - for i from 0 <= i < length: - array_data[i] = func(state, (<double *>(itera.dataptr))[0]) - PyArray_ITER_NEXT(itera) + with lock, nogil: + for i from 0 <= i < length: + array_data[i] = func(state, (<double *>(itera.dataptr))[0]) + PyArray_ITER_NEXT(itera) else: array = <ndarray>np.empty(size, np.float64) array_data = <double *>PyArray_DATA(array) @@ -184,14 +193,15 @@ cdef object cont1_array(rk_state *state, rk_cont1 func, object size, ndarray oa) <void *>oa) if (multi.size != PyArray_SIZE(array)): raise ValueError("size is not compatible with inputs") - for i from 0 <= i < multi.size: - oa_data = <double *>PyArray_MultiIter_DATA(multi, 1) - array_data[i] = func(state, oa_data[0]) - PyArray_MultiIter_NEXTi(multi, 1) + with lock, nogil: + for i from 0 <= i < multi.size: + oa_data = <double *>PyArray_MultiIter_DATA(multi, 1) + array_data[i] = func(state, oa_data[0]) + PyArray_MultiIter_NEXTi(multi, 1) return array cdef object cont2_array_sc(rk_state *state, rk_cont2 func, object size, double a, - double b): + double b, object lock): cdef double *array_data cdef ndarray array "arrayObject" cdef npy_intp length @@ -203,13 +213,14 @@ cdef object cont2_array_sc(rk_state *state, rk_cont2 func, object size, double a array = <ndarray>np.empty(size, np.float64) length = PyArray_SIZE(array) array_data = <double *>PyArray_DATA(array) - for i from 0 <= i < length: - array_data[i] = func(state, a, b) + with lock, nogil: + for i from 0 <= i < length: + array_data[i] = func(state, a, b) return array cdef object cont2_array(rk_state *state, rk_cont2 func, object size, - ndarray oa, ndarray ob): + ndarray oa, ndarray ob, object lock): cdef double *array_data cdef double *oa_data cdef double *ob_data @@ -222,27 +233,29 @@ cdef object cont2_array(rk_state *state, rk_cont2 func, object size, multi = <broadcast> PyArray_MultiIterNew(2, <void *>oa, <void *>ob) array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_DOUBLE) array_data = <double *>PyArray_DATA(array) - for i from 0 <= i < multi.size: - oa_data = <double *>PyArray_MultiIter_DATA(multi, 0) - ob_data = <double *>PyArray_MultiIter_DATA(multi, 1) - array_data[i] = func(state, oa_data[0], ob_data[0]) - PyArray_MultiIter_NEXT(multi) + with lock, nogil: + for i from 0 <= i < multi.size: + oa_data = <double *>PyArray_MultiIter_DATA(multi, 0) + ob_data = <double *>PyArray_MultiIter_DATA(multi, 1) + array_data[i] = func(state, oa_data[0], ob_data[0]) + PyArray_MultiIter_NEXT(multi) else: array = <ndarray>np.empty(size, np.float64) array_data = <double *>PyArray_DATA(array) multi = <broadcast>PyArray_MultiIterNew(3, <void*>array, <void *>oa, <void *>ob) if (multi.size != PyArray_SIZE(array)): raise ValueError("size is not compatible with inputs") - for i from 0 <= i < multi.size: - oa_data = <double *>PyArray_MultiIter_DATA(multi, 1) - ob_data = <double *>PyArray_MultiIter_DATA(multi, 2) - array_data[i] = func(state, oa_data[0], ob_data[0]) - PyArray_MultiIter_NEXTi(multi, 1) - PyArray_MultiIter_NEXTi(multi, 2) + with lock, nogil: + for i from 0 <= i < multi.size: + oa_data = <double *>PyArray_MultiIter_DATA(multi, 1) + ob_data = <double *>PyArray_MultiIter_DATA(multi, 2) + array_data[i] = func(state, oa_data[0], ob_data[0]) + PyArray_MultiIter_NEXTi(multi, 1) + PyArray_MultiIter_NEXTi(multi, 2) return array cdef object cont3_array_sc(rk_state *state, rk_cont3 func, object size, double a, - double b, double c): + double b, double c, object lock): cdef double *array_data cdef ndarray array "arrayObject" @@ -255,12 +268,13 @@ cdef object cont3_array_sc(rk_state *state, rk_cont3 func, object size, double a array = <ndarray>np.empty(size, np.float64) length = PyArray_SIZE(array) array_data = <double *>PyArray_DATA(array) - for i from 0 <= i < length: - array_data[i] = func(state, a, b, c) + with lock, nogil: + for i from 0 <= i < length: + array_data[i] = func(state, a, b, c) return array -cdef object cont3_array(rk_state *state, rk_cont3 func, object size, ndarray oa, - ndarray ob, ndarray oc): +cdef object cont3_array(rk_state *state, rk_cont3 func, object size, + ndarray oa, ndarray ob, ndarray oc, object lock): cdef double *array_data cdef double *oa_data @@ -275,12 +289,13 @@ cdef object cont3_array(rk_state *state, rk_cont3 func, object size, ndarray oa, multi = <broadcast> PyArray_MultiIterNew(3, <void *>oa, <void *>ob, <void *>oc) array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_DOUBLE) array_data = <double *>PyArray_DATA(array) - for i from 0 <= i < multi.size: - oa_data = <double *>PyArray_MultiIter_DATA(multi, 0) - ob_data = <double *>PyArray_MultiIter_DATA(multi, 1) - oc_data = <double *>PyArray_MultiIter_DATA(multi, 2) - array_data[i] = func(state, oa_data[0], ob_data[0], oc_data[0]) - PyArray_MultiIter_NEXT(multi) + with lock, nogil: + for i from 0 <= i < multi.size: + oa_data = <double *>PyArray_MultiIter_DATA(multi, 0) + ob_data = <double *>PyArray_MultiIter_DATA(multi, 1) + oc_data = <double *>PyArray_MultiIter_DATA(multi, 2) + array_data[i] = func(state, oa_data[0], ob_data[0], oc_data[0]) + PyArray_MultiIter_NEXT(multi) else: array = <ndarray>np.empty(size, np.float64) array_data = <double *>PyArray_DATA(array) @@ -288,15 +303,16 @@ cdef object cont3_array(rk_state *state, rk_cont3 func, object size, ndarray oa, <void *>ob, <void *>oc) if (multi.size != PyArray_SIZE(array)): raise ValueError("size is not compatible with inputs") - for i from 0 <= i < multi.size: - oa_data = <double *>PyArray_MultiIter_DATA(multi, 1) - ob_data = <double *>PyArray_MultiIter_DATA(multi, 2) - oc_data = <double *>PyArray_MultiIter_DATA(multi, 3) - array_data[i] = func(state, oa_data[0], ob_data[0], oc_data[0]) - PyArray_MultiIter_NEXT(multi) + with lock, nogil: + for i from 0 <= i < multi.size: + oa_data = <double *>PyArray_MultiIter_DATA(multi, 1) + ob_data = <double *>PyArray_MultiIter_DATA(multi, 2) + oc_data = <double *>PyArray_MultiIter_DATA(multi, 3) + array_data[i] = func(state, oa_data[0], ob_data[0], oc_data[0]) + PyArray_MultiIter_NEXT(multi) return array -cdef object disc0_array(rk_state *state, rk_disc0 func, object size): +cdef object disc0_array(rk_state *state, rk_disc0 func, object size, object lock): cdef long *array_data cdef ndarray array "arrayObject" cdef npy_intp length @@ -308,11 +324,13 @@ cdef object disc0_array(rk_state *state, rk_disc0 func, object size): array = <ndarray>np.empty(size, int) length = PyArray_SIZE(array) array_data = <long *>PyArray_DATA(array) - for i from 0 <= i < length: - array_data[i] = func(state) + with lock, nogil: + for i from 0 <= i < length: + array_data[i] = func(state) return array -cdef object discnp_array_sc(rk_state *state, rk_discnp func, object size, long n, double p): +cdef object discnp_array_sc(rk_state *state, rk_discnp func, object size, + long n, double p, object lock): cdef long *array_data cdef ndarray array "arrayObject" cdef npy_intp length @@ -324,11 +342,13 @@ cdef object discnp_array_sc(rk_state *state, rk_discnp func, object size, long n array = <ndarray>np.empty(size, int) length = PyArray_SIZE(array) array_data = <long *>PyArray_DATA(array) - for i from 0 <= i < length: - array_data[i] = func(state, n, p) + with lock, nogil: + for i from 0 <= i < length: + array_data[i] = func(state, n, p) return array -cdef object discnp_array(rk_state *state, rk_discnp func, object size, ndarray on, ndarray op): +cdef object discnp_array(rk_state *state, rk_discnp func, object size, + ndarray on, ndarray op, object lock): cdef long *array_data cdef ndarray array "arrayObject" cdef npy_intp length @@ -341,27 +361,30 @@ cdef object discnp_array(rk_state *state, rk_discnp func, object size, ndarray o multi = <broadcast> PyArray_MultiIterNew(2, <void *>on, <void *>op) array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_LONG) array_data = <long *>PyArray_DATA(array) - for i from 0 <= i < multi.size: - on_data = <long *>PyArray_MultiIter_DATA(multi, 0) - op_data = <double *>PyArray_MultiIter_DATA(multi, 1) - array_data[i] = func(state, on_data[0], op_data[0]) - PyArray_MultiIter_NEXT(multi) + with lock, nogil: + for i from 0 <= i < multi.size: + on_data = <long *>PyArray_MultiIter_DATA(multi, 0) + op_data = <double *>PyArray_MultiIter_DATA(multi, 1) + array_data[i] = func(state, on_data[0], op_data[0]) + PyArray_MultiIter_NEXT(multi) else: array = <ndarray>np.empty(size, int) array_data = <long *>PyArray_DATA(array) multi = <broadcast>PyArray_MultiIterNew(3, <void*>array, <void *>on, <void *>op) if (multi.size != PyArray_SIZE(array)): raise ValueError("size is not compatible with inputs") - for i from 0 <= i < multi.size: - on_data = <long *>PyArray_MultiIter_DATA(multi, 1) - op_data = <double *>PyArray_MultiIter_DATA(multi, 2) - array_data[i] = func(state, on_data[0], op_data[0]) - PyArray_MultiIter_NEXTi(multi, 1) - PyArray_MultiIter_NEXTi(multi, 2) + with lock, nogil: + for i from 0 <= i < multi.size: + on_data = <long *>PyArray_MultiIter_DATA(multi, 1) + op_data = <double *>PyArray_MultiIter_DATA(multi, 2) + array_data[i] = func(state, on_data[0], op_data[0]) + PyArray_MultiIter_NEXTi(multi, 1) + PyArray_MultiIter_NEXTi(multi, 2) return array -cdef object discdd_array_sc(rk_state *state, rk_discdd func, object size, double n, double p): +cdef object discdd_array_sc(rk_state *state, rk_discdd func, object size, + double n, double p, object lock): cdef long *array_data cdef ndarray array "arrayObject" cdef npy_intp length @@ -373,11 +396,13 @@ cdef object discdd_array_sc(rk_state *state, rk_discdd func, object size, double array = <ndarray>np.empty(size, int) length = PyArray_SIZE(array) array_data = <long *>PyArray_DATA(array) - for i from 0 <= i < length: - array_data[i] = func(state, n, p) + with lock, nogil: + for i from 0 <= i < length: + array_data[i] = func(state, n, p) return array -cdef object discdd_array(rk_state *state, rk_discdd func, object size, ndarray on, ndarray op): +cdef object discdd_array(rk_state *state, rk_discdd func, object size, + ndarray on, ndarray op, object lock): cdef long *array_data cdef ndarray array "arrayObject" cdef npy_intp length @@ -390,28 +415,30 @@ cdef object discdd_array(rk_state *state, rk_discdd func, object size, ndarray o multi = <broadcast> PyArray_MultiIterNew(2, <void *>on, <void *>op) array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_LONG) array_data = <long *>PyArray_DATA(array) - for i from 0 <= i < multi.size: - on_data = <double *>PyArray_MultiIter_DATA(multi, 0) - op_data = <double *>PyArray_MultiIter_DATA(multi, 1) - array_data[i] = func(state, on_data[0], op_data[0]) - PyArray_MultiIter_NEXT(multi) + with lock, nogil: + for i from 0 <= i < multi.size: + on_data = <double *>PyArray_MultiIter_DATA(multi, 0) + op_data = <double *>PyArray_MultiIter_DATA(multi, 1) + array_data[i] = func(state, on_data[0], op_data[0]) + PyArray_MultiIter_NEXT(multi) else: array = <ndarray>np.empty(size, int) array_data = <long *>PyArray_DATA(array) multi = <broadcast>PyArray_MultiIterNew(3, <void*>array, <void *>on, <void *>op) if (multi.size != PyArray_SIZE(array)): raise ValueError("size is not compatible with inputs") - for i from 0 <= i < multi.size: - on_data = <double *>PyArray_MultiIter_DATA(multi, 1) - op_data = <double *>PyArray_MultiIter_DATA(multi, 2) - array_data[i] = func(state, on_data[0], op_data[0]) - PyArray_MultiIter_NEXTi(multi, 1) - PyArray_MultiIter_NEXTi(multi, 2) + with lock, nogil: + for i from 0 <= i < multi.size: + on_data = <double *>PyArray_MultiIter_DATA(multi, 1) + op_data = <double *>PyArray_MultiIter_DATA(multi, 2) + array_data[i] = func(state, on_data[0], op_data[0]) + PyArray_MultiIter_NEXTi(multi, 1) + PyArray_MultiIter_NEXTi(multi, 2) return array cdef object discnmN_array_sc(rk_state *state, rk_discnmN func, object size, - long n, long m, long N): + long n, long m, long N, object lock): cdef long *array_data cdef ndarray array "arrayObject" cdef npy_intp length @@ -423,12 +450,13 @@ cdef object discnmN_array_sc(rk_state *state, rk_discnmN func, object size, array = <ndarray>np.empty(size, int) length = PyArray_SIZE(array) array_data = <long *>PyArray_DATA(array) - for i from 0 <= i < length: - array_data[i] = func(state, n, m, N) + with lock, nogil: + for i from 0 <= i < length: + array_data[i] = func(state, n, m, N) return array cdef object discnmN_array(rk_state *state, rk_discnmN func, object size, - ndarray on, ndarray om, ndarray oN): + ndarray on, ndarray om, ndarray oN, object lock): cdef long *array_data cdef long *on_data cdef long *om_data @@ -442,12 +470,13 @@ cdef object discnmN_array(rk_state *state, rk_discnmN func, object size, multi = <broadcast> PyArray_MultiIterNew(3, <void *>on, <void *>om, <void *>oN) array = <ndarray> PyArray_SimpleNew(multi.nd, multi.dimensions, NPY_LONG) array_data = <long *>PyArray_DATA(array) - for i from 0 <= i < multi.size: - on_data = <long *>PyArray_MultiIter_DATA(multi, 0) - om_data = <long *>PyArray_MultiIter_DATA(multi, 1) - oN_data = <long *>PyArray_MultiIter_DATA(multi, 2) - array_data[i] = func(state, on_data[0], om_data[0], oN_data[0]) - PyArray_MultiIter_NEXT(multi) + with lock, nogil: + for i from 0 <= i < multi.size: + on_data = <long *>PyArray_MultiIter_DATA(multi, 0) + om_data = <long *>PyArray_MultiIter_DATA(multi, 1) + oN_data = <long *>PyArray_MultiIter_DATA(multi, 2) + array_data[i] = func(state, on_data[0], om_data[0], oN_data[0]) + PyArray_MultiIter_NEXT(multi) else: array = <ndarray>np.empty(size, int) array_data = <long *>PyArray_DATA(array) @@ -455,16 +484,18 @@ cdef object discnmN_array(rk_state *state, rk_discnmN func, object size, <void *>oN) if (multi.size != PyArray_SIZE(array)): raise ValueError("size is not compatible with inputs") - for i from 0 <= i < multi.size: - on_data = <long *>PyArray_MultiIter_DATA(multi, 1) - om_data = <long *>PyArray_MultiIter_DATA(multi, 2) - oN_data = <long *>PyArray_MultiIter_DATA(multi, 3) - array_data[i] = func(state, on_data[0], om_data[0], oN_data[0]) - PyArray_MultiIter_NEXT(multi) + with lock, nogil: + for i from 0 <= i < multi.size: + on_data = <long *>PyArray_MultiIter_DATA(multi, 1) + om_data = <long *>PyArray_MultiIter_DATA(multi, 2) + oN_data = <long *>PyArray_MultiIter_DATA(multi, 3) + array_data[i] = func(state, on_data[0], om_data[0], oN_data[0]) + PyArray_MultiIter_NEXT(multi) return array -cdef object discd_array_sc(rk_state *state, rk_discd func, object size, double a): +cdef object discd_array_sc(rk_state *state, rk_discd func, object size, + double a, object lock): cdef long *array_data cdef ndarray array "arrayObject" cdef npy_intp length @@ -476,11 +507,13 @@ cdef object discd_array_sc(rk_state *state, rk_discd func, object size, double a array = <ndarray>np.empty(size, int) length = PyArray_SIZE(array) array_data = <long *>PyArray_DATA(array) - for i from 0 <= i < length: - array_data[i] = func(state, a) + with lock, nogil: + for i from 0 <= i < length: + array_data[i] = func(state, a) return array -cdef object discd_array(rk_state *state, rk_discd func, object size, ndarray oa): +cdef object discd_array(rk_state *state, rk_discd func, object size, ndarray oa, + object lock): cdef long *array_data cdef double *oa_data cdef ndarray array "arrayObject" @@ -495,19 +528,21 @@ cdef object discd_array(rk_state *state, rk_discd func, object size, ndarray oa) length = PyArray_SIZE(array) array_data = <long *>PyArray_DATA(array) itera = <flatiter>PyArray_IterNew(<object>oa) - for i from 0 <= i < length: - array_data[i] = func(state, (<double *>(itera.dataptr))[0]) - PyArray_ITER_NEXT(itera) + with lock, nogil: + for i from 0 <= i < length: + array_data[i] = func(state, (<double *>(itera.dataptr))[0]) + PyArray_ITER_NEXT(itera) else: array = <ndarray>np.empty(size, int) array_data = <long *>PyArray_DATA(array) multi = <broadcast>PyArray_MultiIterNew(2, <void *>array, <void *>oa) if (multi.size != PyArray_SIZE(array)): raise ValueError("size is not compatible with inputs") - for i from 0 <= i < multi.size: - oa_data = <double *>PyArray_MultiIter_DATA(multi, 1) - array_data[i] = func(state, oa_data[0]) - PyArray_MultiIter_NEXTi(multi, 1) + with lock, nogil: + for i from 0 <= i < multi.size: + oa_data = <double *>PyArray_MultiIter_DATA(multi, 1) + array_data[i] = func(state, oa_data[0]) + PyArray_MultiIter_NEXTi(multi, 1) return array cdef double kahan_sum(double *darr, npy_intp n): @@ -522,6 +557,16 @@ cdef double kahan_sum(double *darr, npy_intp n): sum = t return sum +def _shape_from_size(size, d): + if size is None: + shape = (d,) + else: + try: + shape = (operator.index(size), d) + except TypeError: + shape = tuple(size) + (d,) + return shape + cdef class RandomState: """ RandomState(seed=None) @@ -556,12 +601,14 @@ cdef class RandomState: """ cdef rk_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): self.internal_state = <rk_state*>PyMem_Malloc(sizeof(rk_state)) self.seed(seed) + self.lock = Lock() def __dealloc__(self): if self.internal_state != NULL: @@ -581,6 +628,7 @@ cdef class RandomState: ---------- seed : int or array_like, optional Seed for `RandomState`. + Must be convertable to 32 bit unsigned integers. See Also -------- @@ -589,15 +637,19 @@ cdef class RandomState: """ cdef rk_error errcode cdef ndarray obj "arrayObject_obj" - if seed is None: - errcode = rk_randomseed(self.internal_state) - elif type(seed) is int: - rk_seed(seed, self.internal_state) - elif isinstance(seed, np.integer): - iseed = int(seed) - rk_seed(iseed, self.internal_state) - else: - obj = <ndarray>PyArray_ContiguousFromObject(seed, NPY_LONG, 1, 1) + try: + if seed is None: + errcode = rk_randomseed(self.internal_state) + else: + idx = operator.index(seed) + if idx > int(2**32 - 1) or idx < 0: + raise ValueError("Seed must be between 0 and 4294967295") + rk_seed(idx, self.internal_state) + except TypeError: + obj = np.asarray(seed).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('L', casting='unsafe') init_by_array(self.internal_state, <unsigned long *>PyArray_DATA(obj), PyArray_DIM(obj, 0)) @@ -734,8 +786,9 @@ cdef class RandomState: Parameters ---------- size : int or tuple of ints, optional - Defines the shape of the returned array of random floats. If None - (the default), returns a single float. + 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 ------- @@ -760,7 +813,7 @@ cdef class RandomState: [-1.23204345, -1.75224494]]) """ - return cont0_array(self.internal_state, rk_double, size) + return cont0_array(self.internal_state, rk_double, size, self.lock) def tomaxint(self, size=None): """ @@ -773,10 +826,10 @@ cdef class RandomState: Parameters ---------- - size : tuple of ints, int, optional - Shape of output. If this is, for example, (m,n,k), m*n*k samples - are generated. If no shape is specified, a single sample is - returned. + 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 ------- @@ -807,7 +860,7 @@ cdef class RandomState: [ True, True]]], dtype=bool) """ - return disc0_array(self.internal_state, rk_long, size) + return disc0_array(self.internal_state, rk_long, size, self.lock) def randint(self, low, high=None, size=None): """ @@ -829,8 +882,9 @@ cdef class RandomState: If provided, one above the largest (signed) integer to be drawn from the distribution (see above for behavior if ``high=None``). size : int or tuple of ints, optional - Output shape. Default is None, in which case a single int is - returned. + 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 ------- @@ -884,9 +938,10 @@ cdef class RandomState: array = <ndarray>np.empty(size, int) length = PyArray_SIZE(array) array_data = <long *>PyArray_DATA(array) - for i from 0 <= i < length: - rv = lo + <long>rk_interval(diff, self. internal_state) - array_data[i] = rv + with self.lock, nogil: + for i from 0 <= i < length: + rv = lo + <long>rk_interval(diff, self. internal_state) + array_data[i] = rv return array def bytes(self, npy_intp length): @@ -913,7 +968,8 @@ cdef class RandomState: """ cdef void *bytes bytestring = empty_py_bytes(length, &bytes) - rk_fill(bytes, length, self.internal_state) + with self.lock, nogil: + rk_fill(bytes, length, self.internal_state) return bytestring @@ -931,13 +987,14 @@ cdef class RandomState: If an ndarray, a random sample is generated from its elements. If an int, the random sample is generated as if a was np.arange(n) size : int or tuple of ints, optional - Output shape. Default is None, in which case a single value is - returned. + 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. replace : boolean, optional Whether the sample is with or without replacement p : 1-D array-like, optional The probabilities associated with each entry in a. - If not given the sample assumes a uniform distribtion over all + If not given the sample assumes a uniform distribution over all entries in a. Returns @@ -976,7 +1033,7 @@ cdef class RandomState: >>> np.random.choice(5, 3, replace=False) array([3,1,0]) - >>> #This is equivalent to np.random.shuffle(np.arange(5))[:3] + >>> #This is equivalent to np.random.permutation(np.arange(5))[:3] Generate a non-uniform random sample from np.arange(5) of size 3 without replacement: @@ -1012,14 +1069,17 @@ cdef class RandomState: raise ValueError("a must be non-empty") if None != p: - p = np.array(p, dtype=np.double, ndmin=1, copy=False) + d = len(p) + p = <ndarray>PyArray_ContiguousFromObject(p, NPY_DOUBLE, 1, 1) + pix = <double*>PyArray_DATA(p) + if p.ndim != 1: raise ValueError("p must be 1-dimensional") if p.size != pop_size: raise ValueError("a and p must have same size") - if np.any(p < 0): + if np.logical_or.reduce(p < 0): raise ValueError("probabilities are not non-negative") - if not np.allclose(p.sum(), 1): + if abs(kahan_sum(pix, d) - 1.) > 1e-8: raise ValueError("probabilities do not sum to 1") shape = size @@ -1044,7 +1104,7 @@ cdef class RandomState: "population when 'replace=False'") if None != p: - if np.sum(p > 0) < size: + if np.count_nonzero(p > 0) < size: raise ValueError("Fewer non-zero entries in p than size") n_uniq = 0 p = p.copy() @@ -1091,7 +1151,7 @@ cdef class RandomState: def uniform(self, low=0.0, high=1.0, size=None): """ - uniform(low=0.0, high=1.0, size=1) + uniform(low=0.0, high=1.0, size=None) Draw samples from a uniform distribution. @@ -1109,9 +1169,9 @@ cdef class RandomState: Upper boundary of the output interval. All values generated will be less than high. The default value is 1.0. size : int or tuple of ints, optional - Shape of output. If the given size is, for example, (m,n,k), - m*n*k samples are generated. If no shape is specified, a single sample - is returned. + 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 ------- @@ -1166,7 +1226,9 @@ cdef class RandomState: flow = PyFloat_AsDouble(low) fhigh = PyFloat_AsDouble(high) if not PyErr_Occurred(): - return cont2_array_sc(self.internal_state, rk_uniform, size, flow, fhigh-flow) + return cont2_array_sc(self.internal_state, rk_uniform, size, flow, + fhigh-flow, self.lock) + PyErr_Clear() olow = <ndarray>PyArray_FROM_OTF(low, NPY_DOUBLE, NPY_ARRAY_ALIGNED) ohigh = <ndarray>PyArray_FROM_OTF(high, NPY_DOUBLE, NPY_ARRAY_ALIGNED) @@ -1174,7 +1236,8 @@ cdef class RandomState: Py_INCREF(temp) # needed to get around Pyrex's automatic reference-counting # rules because EnsureArray steals a reference odiff = <ndarray>PyArray_EnsureArray(temp) - return cont2_array(self.internal_state, rk_uniform, size, olow, odiff) + return cont2_array(self.internal_state, rk_uniform, size, olow, odiff, + self.lock) def rand(self, *args): """ @@ -1297,7 +1360,9 @@ cdef class RandomState: If provided, the largest (signed) integer to be drawn from the distribution (see above for behavior if ``high=None``). size : int or tuple of ints, optional - Output shape. Default is None, in which case a single int is returned. + 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 ------- @@ -1364,8 +1429,9 @@ cdef class RandomState: Parameters ---------- size : int or tuple of ints, optional - Output shape. Default is None, in which case a single value is - returned. + 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 ------- @@ -1385,7 +1451,7 @@ cdef class RandomState: (3, 4, 2) """ - return cont0_array(self.internal_state, rk_gauss, size) + return cont0_array(self.internal_state, rk_gauss, size, self.lock) def normal(self, loc=0.0, scale=1.0, size=None): """ @@ -1409,9 +1475,10 @@ cdef class RandomState: Mean ("centre") of the distribution. scale : float Standard deviation (spread or "width") of the distribution. - size : tuple of ints + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. See Also -------- @@ -1477,7 +1544,8 @@ cdef class RandomState: if not PyErr_Occurred(): if fscale <= 0: raise ValueError("scale <= 0") - return cont2_array_sc(self.internal_state, rk_normal, size, floc, fscale) + return cont2_array_sc(self.internal_state, rk_normal, size, floc, + fscale, self.lock) PyErr_Clear() @@ -1485,7 +1553,8 @@ cdef class RandomState: oscale = <ndarray>PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(oscale, 0)): raise ValueError("scale <= 0") - return cont2_array(self.internal_state, rk_normal, size, oloc, oscale) + return cont2_array(self.internal_state, rk_normal, size, oloc, oscale, + self.lock) def beta(self, a, b, size=None): """ @@ -1513,9 +1582,10 @@ cdef class RandomState: Alpha, non-negative. b : float Beta, non-negative. - size : tuple of ints, optional - The number of samples to draw. The output is packed according to - the size given. + 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 ------- @@ -1534,7 +1604,8 @@ cdef class RandomState: raise ValueError("a <= 0") if fb <= 0: raise ValueError("b <= 0") - return cont2_array_sc(self.internal_state, rk_beta, size, fa, fb) + return cont2_array_sc(self.internal_state, rk_beta, size, fa, fb, + self.lock) PyErr_Clear() @@ -1544,7 +1615,8 @@ cdef class RandomState: raise ValueError("a <= 0") if np.any(np.less_equal(ob, 0)): raise ValueError("b <= 0") - return cont2_array(self.internal_state, rk_beta, size, oa, ob) + return cont2_array(self.internal_state, rk_beta, size, oa, ob, + self.lock) def exponential(self, scale=1.0, size=None): """ @@ -1570,9 +1642,10 @@ cdef class RandomState: ---------- scale : float The scale parameter, :math:`\\beta = 1/\\lambda`. - size : tuple of ints - Number of samples to draw. The output is shaped - according to `size`. + 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. References ---------- @@ -1591,14 +1664,16 @@ cdef class RandomState: if not PyErr_Occurred(): if fscale <= 0: raise ValueError("scale <= 0") - return cont1_array_sc(self.internal_state, rk_exponential, size, fscale) + return cont1_array_sc(self.internal_state, rk_exponential, size, + fscale, self.lock) PyErr_Clear() oscale = <ndarray> PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(oscale, 0.0)): raise ValueError("scale <= 0") - return cont1_array(self.internal_state, rk_exponential, size, oscale) + return cont1_array(self.internal_state, rk_exponential, size, oscale, + self.lock) def standard_exponential(self, size=None): """ @@ -1611,8 +1686,10 @@ cdef class RandomState: Parameters ---------- - size : int or tuple of ints - Shape of the output. + 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 ------- @@ -1626,7 +1703,8 @@ cdef class RandomState: >>> n = np.random.standard_exponential((3, 8000)) """ - return cont0_array(self.internal_state, rk_standard_exponential, size) + return cont0_array(self.internal_state, rk_standard_exponential, size, + self.lock) def standard_gamma(self, shape, size=None): """ @@ -1641,9 +1719,10 @@ cdef class RandomState: ---------- shape : float Parameter, should be > 0. - size : int or tuple of ints + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. Returns ------- @@ -1702,13 +1781,14 @@ cdef class RandomState: if not PyErr_Occurred(): if fshape <= 0: raise ValueError("shape <= 0") - return cont1_array_sc(self.internal_state, rk_standard_gamma, size, fshape) + return cont1_array_sc(self.internal_state, rk_standard_gamma, size, fshape, self.lock) PyErr_Clear() oshape = <ndarray> PyArray_FROM_OTF(shape, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(oshape, 0.0)): raise ValueError("shape <= 0") - return cont1_array(self.internal_state, rk_standard_gamma, size, oshape) + return cont1_array(self.internal_state, rk_standard_gamma, size, + oshape, self.lock) def gamma(self, shape, scale=1.0, size=None): """ @@ -1726,9 +1806,10 @@ cdef class RandomState: The shape of the gamma distribution. scale : scalar > 0, optional The scale of the gamma distribution. Default is equal to 1. - size : shape_tuple, optional + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. Returns ------- @@ -1790,7 +1871,8 @@ cdef class RandomState: raise ValueError("shape <= 0") if fscale <= 0: raise ValueError("scale <= 0") - return cont2_array_sc(self.internal_state, rk_gamma, size, fshape, fscale) + return cont2_array_sc(self.internal_state, rk_gamma, size, fshape, + fscale, self.lock) PyErr_Clear() oshape = <ndarray>PyArray_FROM_OTF(shape, NPY_DOUBLE, NPY_ARRAY_ALIGNED) @@ -1799,7 +1881,8 @@ cdef class RandomState: raise ValueError("shape <= 0") if np.any(np.less_equal(oscale, 0.0)): raise ValueError("scale <= 0") - return cont2_array(self.internal_state, rk_gamma, size, oshape, oscale) + return cont2_array(self.internal_state, rk_gamma, size, oshape, oscale, + self.lock) def f(self, dfnum, dfden, size=None): """ @@ -1822,10 +1905,10 @@ cdef class RandomState: Degrees of freedom in numerator. Should be greater than zero. dfden : float Degrees of freedom in denominator. Should be greater than zero. - size : {tuple, int}, optional - Output shape. If the given shape is, e.g., ``(m, n, k)``, - then ``m * n * k`` samples are drawn. By default only one sample - is returned. + 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 ------- @@ -1891,7 +1974,8 @@ cdef class RandomState: raise ValueError("shape <= 0") if fdfden <= 0: raise ValueError("scale <= 0") - return cont2_array_sc(self.internal_state, rk_f, size, fdfnum, fdfden) + return cont2_array_sc(self.internal_state, rk_f, size, fdfnum, + fdfden, self.lock) PyErr_Clear() @@ -1901,7 +1985,8 @@ cdef class RandomState: raise ValueError("dfnum <= 0") if np.any(np.less_equal(odfden, 0.0)): raise ValueError("dfden <= 0") - return cont2_array(self.internal_state, rk_f, size, odfnum, odfden) + return cont2_array(self.internal_state, rk_f, size, odfnum, odfden, + self.lock) def noncentral_f(self, dfnum, dfden, nonc, size=None): """ @@ -1922,9 +2007,10 @@ cdef class RandomState: Parameter, should be > 1. nonc : float Parameter, should be >= 0. - size : int or tuple of ints - Output shape. If the given shape is, e.g., ``(m, n, k)``, then - ``m * n * k`` samples are drawn. + 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 ------- @@ -1981,7 +2067,7 @@ cdef class RandomState: if fnonc < 0: raise ValueError("nonc < 0") return cont3_array_sc(self.internal_state, rk_noncentral_f, size, - fdfnum, fdfden, fnonc) + fdfnum, fdfden, fnonc, self.lock) PyErr_Clear() @@ -1996,7 +2082,7 @@ cdef class RandomState: if np.any(np.less(ononc, 0.0)): raise ValueError("nonc < 0") return cont3_array(self.internal_state, rk_noncentral_f, size, odfnum, - odfden, ononc) + odfden, ononc, self.lock) def chisquare(self, df, size=None): """ @@ -2013,9 +2099,10 @@ cdef class RandomState: ---------- df : int Number of degrees of freedom. - size : tuple of ints, int, optional - Size of the returned array. By default, a scalar is - returned. + 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 ------- @@ -2067,14 +2154,16 @@ cdef class RandomState: if not PyErr_Occurred(): if fdf <= 0: raise ValueError("df <= 0") - return cont1_array_sc(self.internal_state, rk_chisquare, size, fdf) + return cont1_array_sc(self.internal_state, rk_chisquare, size, fdf, + self.lock) PyErr_Clear() odf = <ndarray>PyArray_FROM_OTF(df, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(odf, 0.0)): raise ValueError("df <= 0") - return cont1_array(self.internal_state, rk_chisquare, size, odf) + return cont1_array(self.internal_state, rk_chisquare, size, odf, + self.lock) def noncentral_chisquare(self, df, nonc, size=None): """ @@ -2091,8 +2180,10 @@ cdef class RandomState: Degrees of freedom, should be >= 1. nonc : float Non-centrality, should be > 0. - size : int or tuple of ints - Shape of the output. + 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. Notes ----- @@ -2155,7 +2246,7 @@ cdef class RandomState: if fnonc <= 0: raise ValueError("nonc <= 0") return cont2_array_sc(self.internal_state, rk_noncentral_chisquare, - size, fdf, fnonc) + size, fdf, fnonc, self.lock) PyErr_Clear() @@ -2166,7 +2257,7 @@ cdef class RandomState: if np.any(np.less_equal(ononc, 0.0)): raise ValueError("nonc < 0") return cont2_array(self.internal_state, rk_noncentral_chisquare, size, - odf, ononc) + odf, ononc, self.lock) def standard_cauchy(self, size=None): """ @@ -2178,8 +2269,10 @@ cdef class RandomState: Parameters ---------- - size : int or tuple of ints - Shape of the output. + 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 ------- @@ -2227,7 +2320,8 @@ cdef class RandomState: >>> plt.show() """ - return cont0_array(self.internal_state, rk_standard_cauchy, size) + return cont0_array(self.internal_state, rk_standard_cauchy, size, + self.lock) def standard_t(self, df, size=None): """ @@ -2244,8 +2338,9 @@ cdef class RandomState: df : int Degrees of freedom, should be > 0. size : int or tuple of ints, optional - Output shape. Default is None, in which case a single value is - returned. + 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 ------- @@ -2321,14 +2416,16 @@ cdef class RandomState: if not PyErr_Occurred(): if fdf <= 0: raise ValueError("df <= 0") - return cont1_array_sc(self.internal_state, rk_standard_t, size, fdf) + return cont1_array_sc(self.internal_state, rk_standard_t, size, + fdf, self.lock) PyErr_Clear() odf = <ndarray> PyArray_FROM_OTF(df, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(odf, 0.0)): raise ValueError("df <= 0") - return cont1_array(self.internal_state, rk_standard_t, size, odf) + return cont1_array(self.internal_state, rk_standard_t, size, odf, + self.lock) def vonmises(self, mu, kappa, size=None): """ @@ -2350,9 +2447,10 @@ cdef class RandomState: Mode ("center") of the distribution. kappa : float Dispersion of the distribution, has to be >=0. - size : int or tuple of int + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. Returns ------- @@ -2414,7 +2512,8 @@ cdef class RandomState: if not PyErr_Occurred(): if fkappa < 0: raise ValueError("kappa < 0") - return cont2_array_sc(self.internal_state, rk_vonmises, size, fmu, fkappa) + return cont2_array_sc(self.internal_state, rk_vonmises, size, fmu, + fkappa, self.lock) PyErr_Clear() @@ -2422,7 +2521,8 @@ cdef class RandomState: okappa = <ndarray> PyArray_FROM_OTF(kappa, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less(okappa, 0.0)): raise ValueError("kappa < 0") - return cont2_array(self.internal_state, rk_vonmises, size, omu, okappa) + return cont2_array(self.internal_state, rk_vonmises, size, omu, okappa, + self.lock) def pareto(self, a, size=None): """ @@ -2448,9 +2548,10 @@ cdef class RandomState: ---------- shape : float, > 0. Shape of the distribution. - size : tuple of ints + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. See Also -------- @@ -2511,14 +2612,15 @@ cdef class RandomState: if not PyErr_Occurred(): if fa <= 0: raise ValueError("a <= 0") - return cont1_array_sc(self.internal_state, rk_pareto, size, fa) + return cont1_array_sc(self.internal_state, rk_pareto, size, fa, + self.lock) PyErr_Clear() oa = <ndarray>PyArray_FROM_OTF(a, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(oa, 0.0)): raise ValueError("a <= 0") - return cont1_array(self.internal_state, rk_pareto, size, oa) + return cont1_array(self.internal_state, rk_pareto, size, oa, self.lock) def weibull(self, a, size=None): """ @@ -2540,9 +2642,10 @@ cdef class RandomState: ---------- a : float Shape of the distribution. - size : tuple of ints + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. See Also -------- @@ -2611,14 +2714,16 @@ cdef class RandomState: if not PyErr_Occurred(): if fa <= 0: raise ValueError("a <= 0") - return cont1_array_sc(self.internal_state, rk_weibull, size, fa) + return cont1_array_sc(self.internal_state, rk_weibull, size, fa, + self.lock) PyErr_Clear() oa = <ndarray>PyArray_FROM_OTF(a, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(oa, 0.0)): raise ValueError("a <= 0") - return cont1_array(self.internal_state, rk_weibull, size, oa) + return cont1_array(self.internal_state, rk_weibull, size, oa, + self.lock) def power(self, a, size=None): """ @@ -2633,9 +2738,10 @@ cdef class RandomState: ---------- a : float parameter, > 0 - size : tuple of ints + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. Returns ------- @@ -2720,14 +2826,15 @@ cdef class RandomState: if not PyErr_Occurred(): if fa <= 0: raise ValueError("a <= 0") - return cont1_array_sc(self.internal_state, rk_power, size, fa) + return cont1_array_sc(self.internal_state, rk_power, size, fa, + self.lock) PyErr_Clear() oa = <ndarray>PyArray_FROM_OTF(a, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(oa, 0.0)): raise ValueError("a <= 0") - return cont1_array(self.internal_state, rk_power, size, oa) + return cont1_array(self.internal_state, rk_power, size, oa, self.lock) def laplace(self, loc=0.0, scale=1.0, size=None): """ @@ -2747,6 +2854,10 @@ cdef class RandomState: The position, :math:`\\mu`, of the distribution peak. scale : float :math:`\\lambda`, the exponential decay. + 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. Notes ----- @@ -2810,14 +2921,16 @@ cdef class RandomState: if not PyErr_Occurred(): if fscale <= 0: raise ValueError("scale <= 0") - return cont2_array_sc(self.internal_state, rk_laplace, size, floc, fscale) + return cont2_array_sc(self.internal_state, rk_laplace, size, floc, + fscale, self.lock) PyErr_Clear() oloc = PyArray_FROM_OTF(loc, NPY_DOUBLE, NPY_ARRAY_ALIGNED) oscale = PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(oscale, 0.0)): raise ValueError("scale <= 0") - return cont2_array(self.internal_state, rk_laplace, size, oloc, oscale) + return cont2_array(self.internal_state, rk_laplace, size, oloc, oscale, + self.lock) def gumbel(self, loc=0.0, scale=1.0, size=None): """ @@ -2835,9 +2948,10 @@ cdef class RandomState: The location of the mode of the distribution. scale : float The scale parameter of the distribution. - size : tuple of ints + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. Returns ------- @@ -2941,14 +3055,16 @@ cdef class RandomState: if not PyErr_Occurred(): if fscale <= 0: raise ValueError("scale <= 0") - return cont2_array_sc(self.internal_state, rk_gumbel, size, floc, fscale) + return cont2_array_sc(self.internal_state, rk_gumbel, size, floc, + fscale, self.lock) PyErr_Clear() oloc = PyArray_FROM_OTF(loc, NPY_DOUBLE, NPY_ARRAY_ALIGNED) oscale = PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(oscale, 0.0)): raise ValueError("scale <= 0") - return cont2_array(self.internal_state, rk_gumbel, size, oloc, oscale) + return cont2_array(self.internal_state, rk_gumbel, size, oloc, oscale, + self.lock) def logistic(self, loc=0.0, scale=1.0, size=None): """ @@ -2965,9 +3081,10 @@ cdef class RandomState: scale : float > 0. - size : {tuple, int} + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. Returns ------- @@ -3029,14 +3146,16 @@ cdef class RandomState: if not PyErr_Occurred(): if fscale <= 0: raise ValueError("scale <= 0") - return cont2_array_sc(self.internal_state, rk_logistic, size, floc, fscale) + return cont2_array_sc(self.internal_state, rk_logistic, size, floc, + fscale, self.lock) PyErr_Clear() oloc = PyArray_FROM_OTF(loc, NPY_DOUBLE, NPY_ARRAY_ALIGNED) oscale = PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(oscale, 0.0)): raise ValueError("scale <= 0") - return cont2_array(self.internal_state, rk_logistic, size, oloc, oscale) + return cont2_array(self.internal_state, rk_logistic, size, oloc, + oscale, self.lock) def lognormal(self, mean=0.0, sigma=1.0, size=None): """ @@ -3055,9 +3174,10 @@ cdef class RandomState: Mean value of the underlying normal distribution sigma : float, > 0. Standard deviation of the underlying normal distribution - size : tuple of ints + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. Returns ------- @@ -3149,7 +3269,8 @@ cdef class RandomState: if not PyErr_Occurred(): if fsigma <= 0: raise ValueError("sigma <= 0") - return cont2_array_sc(self.internal_state, rk_lognormal, size, fmean, fsigma) + return cont2_array_sc(self.internal_state, rk_lognormal, size, + fmean, fsigma, self.lock) PyErr_Clear() @@ -3157,7 +3278,8 @@ cdef class RandomState: osigma = PyArray_FROM_OTF(sigma, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(osigma, 0.0)): raise ValueError("sigma <= 0.0") - return cont2_array(self.internal_state, rk_lognormal, size, omean, osigma) + return cont2_array(self.internal_state, rk_lognormal, size, omean, + osigma, self.lock) def rayleigh(self, scale=1.0, size=None): """ @@ -3173,8 +3295,9 @@ cdef class RandomState: scale : scalar Scale, also equals the mode. Should be >= 0. size : int or tuple of ints, optional - Shape of the output. Default is None, in which case a single - value is returned. + 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. Notes ----- @@ -3222,14 +3345,16 @@ cdef class RandomState: if not PyErr_Occurred(): if fscale <= 0: raise ValueError("scale <= 0") - return cont1_array_sc(self.internal_state, rk_rayleigh, size, fscale) + return cont1_array_sc(self.internal_state, rk_rayleigh, size, + fscale, self.lock) PyErr_Clear() oscale = <ndarray>PyArray_FROM_OTF(scale, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(oscale, 0.0)): raise ValueError("scale <= 0.0") - return cont1_array(self.internal_state, rk_rayleigh, size, oscale) + return cont1_array(self.internal_state, rk_rayleigh, size, oscale, + self.lock) def wald(self, mean, scale, size=None): """ @@ -3255,8 +3380,9 @@ cdef class RandomState: scale : scalar Scale parameter, should be >= 0. size : int or tuple of ints, optional - Output shape. Default is None, in which case a single value is - returned. + 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 ------- @@ -3304,7 +3430,8 @@ cdef class RandomState: raise ValueError("mean <= 0") if fscale <= 0: raise ValueError("scale <= 0") - return cont2_array_sc(self.internal_state, rk_wald, size, fmean, fscale) + return cont2_array_sc(self.internal_state, rk_wald, size, fmean, + fscale, self.lock) PyErr_Clear() omean = PyArray_FROM_OTF(mean, NPY_DOUBLE, NPY_ARRAY_ALIGNED) @@ -3313,9 +3440,8 @@ cdef class RandomState: raise ValueError("mean <= 0.0") elif np.any(np.less_equal(oscale,0.0)): raise ValueError("scale <= 0.0") - return cont2_array(self.internal_state, rk_wald, size, omean, oscale) - - + return cont2_array(self.internal_state, rk_wald, size, omean, oscale, + self.lock) def triangular(self, left, mode, right, size=None): """ @@ -3337,8 +3463,9 @@ cdef class RandomState: right : scalar Upper limit, should be larger than `left`. size : int or tuple of ints, optional - Output shape. Default is None, in which case a single value is - returned. + 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 ------- @@ -3388,7 +3515,7 @@ cdef class RandomState: if fleft == fright: raise ValueError("left == right") return cont3_array_sc(self.internal_state, rk_triangular, size, fleft, - fmode, fright) + fmode, fright, self.lock) PyErr_Clear() oleft = <ndarray>PyArray_FROM_OTF(left, NPY_DOUBLE, NPY_ARRAY_ALIGNED) @@ -3402,7 +3529,7 @@ cdef class RandomState: if np.any(np.equal(oleft, oright)): raise ValueError("left == right") return cont3_array(self.internal_state, rk_triangular, size, oleft, - omode, oright) + omode, oright, self.lock) # Complicated, discrete distributions: def binomial(self, n, p, size=None): @@ -3422,9 +3549,10 @@ cdef class RandomState: parameter, >= 0. p : float parameter, >= 0 and <=1. - size : {tuple, int} + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. Returns ------- @@ -3499,7 +3627,10 @@ cdef class RandomState: raise ValueError("p < 0") elif fp > 1: raise ValueError("p > 1") - return discnp_array_sc(self.internal_state, rk_binomial, size, ln, fp) + elif np.isnan(fp): + raise ValueError("p is nan") + return discnp_array_sc(self.internal_state, rk_binomial, size, ln, + fp, self.lock) PyErr_Clear() @@ -3511,7 +3642,8 @@ cdef class RandomState: raise ValueError("p < 0") if np.any(np.greater(p, 1)): raise ValueError("p > 1") - return discnp_array(self.internal_state, rk_binomial, size, on, op) + return discnp_array(self.internal_state, rk_binomial, size, on, op, + self.lock) def negative_binomial(self, n, p, size=None): """ @@ -3529,9 +3661,10 @@ cdef class RandomState: Parameter, > 0. p : float Parameter, >= 0 and <=1. - size : int or tuple of ints - Output shape. If the given shape is, e.g., ``(m, n, k)``, then - ``m * n * k`` samples are drawn. + 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 ------- @@ -3593,7 +3726,7 @@ cdef class RandomState: elif fp > 1: raise ValueError("p > 1") return discdd_array_sc(self.internal_state, rk_negative_binomial, - size, fn, fp) + size, fn, fp, self.lock) PyErr_Clear() @@ -3606,7 +3739,7 @@ cdef class RandomState: if np.any(np.greater(p, 1)): raise ValueError("p > 1") return discdd_array(self.internal_state, rk_negative_binomial, size, - on, op) + on, op, self.lock) def poisson(self, lam=1.0, size=None): """ @@ -3619,11 +3752,13 @@ cdef class RandomState: Parameters ---------- - lam : float - Expectation of interval, should be >= 0. + lam : float or sequence of float + Expectation of interval, should be >= 0. A sequence of expectation + intervals must be broadcastable over the requested size. 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. + 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. Notes ----- @@ -3659,6 +3794,10 @@ cdef class RandomState: >>> count, bins, ignored = plt.hist(s, 14, normed=True) >>> plt.show() + Draw each 100 values for lambda 100 and 500: + + >>> s = np.random.poisson(lam=(100., 500.), size=(100, 2)) + """ cdef ndarray olam cdef double flam @@ -3668,7 +3807,8 @@ cdef class RandomState: raise ValueError("lam < 0") if lam > self.poisson_lam_max: raise ValueError("lam value too large") - return discd_array_sc(self.internal_state, rk_poisson, size, flam) + return discd_array_sc(self.internal_state, rk_poisson, size, flam, + self.lock) PyErr_Clear() @@ -3677,7 +3817,7 @@ cdef class RandomState: raise ValueError("lam < 0") if np.any(np.greater(olam, self.poisson_lam_max)): raise ValueError("lam value too large.") - return discd_array(self.internal_state, rk_poisson, size, olam) + return discd_array(self.internal_state, rk_poisson, size, olam, self.lock) def zipf(self, a, size=None): """ @@ -3697,12 +3837,10 @@ cdef class RandomState: ---------- a : float > 1 Distribution parameter. - size : int or tuple of int, optional + 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; a single integer is equivalent in - its result to providing a mono-tuple, i.e., a 1-D array of length - *size* is returned. The default is None, in which case a single - scalar is returned. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. Returns ------- @@ -3758,14 +3896,15 @@ cdef class RandomState: if not PyErr_Occurred(): if fa <= 1.0: raise ValueError("a <= 1.0") - return discd_array_sc(self.internal_state, rk_zipf, size, fa) + return discd_array_sc(self.internal_state, rk_zipf, size, fa, + self.lock) PyErr_Clear() oa = <ndarray>PyArray_FROM_OTF(a, NPY_DOUBLE, NPY_ARRAY_ALIGNED) if np.any(np.less_equal(oa, 1.0)): raise ValueError("a <= 1.0") - return discd_array(self.internal_state, rk_zipf, size, oa) + return discd_array(self.internal_state, rk_zipf, size, oa, self.lock) def geometric(self, p, size=None): """ @@ -3789,9 +3928,10 @@ cdef class RandomState: ---------- p : float The probability of success of an individual trial. - size : tuple of ints - Number of values to draw from the distribution. The output - is shaped according to `size`. + 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 ------- @@ -3821,7 +3961,8 @@ cdef class RandomState: raise ValueError("p < 0.0") if fp > 1.0: raise ValueError("p > 1.0") - return discd_array_sc(self.internal_state, rk_geometric, size, fp) + return discd_array_sc(self.internal_state, rk_geometric, size, fp, + self.lock) PyErr_Clear() @@ -3831,7 +3972,7 @@ cdef class RandomState: raise ValueError("p < 0.0") if np.any(np.greater(op, 1.0)): raise ValueError("p > 1.0") - return discd_array(self.internal_state, rk_geometric, size, op) + return discd_array(self.internal_state, rk_geometric, size, op, self.lock) def hypergeometric(self, ngood, nbad, nsample, size=None): """ @@ -3853,9 +3994,10 @@ cdef class RandomState: nsample : int or array_like Number of items sampled. Must be at least 1 and at most ``ngood + nbad``. - size : int or tuple of int + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. Returns ------- @@ -3933,8 +4075,8 @@ cdef class RandomState: raise ValueError("nsample < 1") if lngood + lnbad < lnsample: raise ValueError("ngood + nbad < nsample") - return discnmN_array_sc(self.internal_state, rk_hypergeometric, size, - lngood, lnbad, lnsample) + return discnmN_array_sc(self.internal_state, rk_hypergeometric, + size, lngood, lnbad, lnsample, self.lock) PyErr_Clear() @@ -3950,7 +4092,7 @@ cdef class RandomState: if np.any(np.less(np.add(ongood, onbad),onsample)): raise ValueError("ngood + nbad < nsample") return discnmN_array(self.internal_state, rk_hypergeometric, size, - ongood, onbad, onsample) + ongood, onbad, onsample, self.lock) def logseries(self, p, size=None): """ @@ -3967,9 +4109,10 @@ cdef class RandomState: scale : float > 0. - size : {tuple, int} + 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. + ``m * n * k`` samples are drawn. Default is None, in which case a + single value is returned. Returns ------- @@ -4035,7 +4178,8 @@ cdef class RandomState: raise ValueError("p <= 0.0") if fp >= 1.0: raise ValueError("p >= 1.0") - return discd_array_sc(self.internal_state, rk_logseries, size, fp) + return discd_array_sc(self.internal_state, rk_logseries, size, fp, + self.lock) PyErr_Clear() @@ -4044,7 +4188,7 @@ cdef class RandomState: raise ValueError("p <= 0.0") if np.any(np.greater_equal(op, 1.0)): raise ValueError("p >= 1.0") - return discd_array(self.internal_state, rk_logseries, size, op) + return discd_array(self.internal_state, rk_logseries, size, op, self.lock) # Multivariate distributions: def multivariate_normal(self, mean, cov, size=None): @@ -4066,7 +4210,7 @@ cdef class RandomState: Mean of the N-dimensional distribution. cov : 2-D array_like, of shape (N, N) Covariance matrix of the distribution. Must be symmetric and - positive semi-definite for "physically meaningful" results. + positive-semidefinite for "physically meaningful" results. size : int or tuple of ints, optional Given a shape of, for example, ``(m,n,k)``, ``m*n*k`` samples are generated, and packed in an `m`-by-`n`-by-`k` arrangement. Because @@ -4138,44 +4282,55 @@ cdef class RandomState: [True, True] """ + from numpy.dual import svd + # Check preconditions on arguments mean = np.array(mean) cov = np.array(cov) if size is None: shape = [] + elif isinstance(size, (int, long, np.integer)): + shape = [size] else: shape = size + if len(mean.shape) != 1: raise ValueError("mean must be 1 dimensional") if (len(cov.shape) != 2) or (cov.shape[0] != cov.shape[1]): raise ValueError("cov must be 2 dimensional and square") if mean.shape[0] != cov.shape[0]: raise ValueError("mean and cov must have same length") - # Compute shape of output - if isinstance(shape, (int, long, np.integer)): - shape = [shape] + + # Compute shape of output and create a matrix of independent + # standard normally distributed random numbers. The matrix has rows + # with the same length as mean and as many rows are necessary to + # form a matrix of shape final_shape. final_shape = list(shape[:]) final_shape.append(mean.shape[0]) - # Create a matrix of independent standard normally distributed random - # numbers. The matrix has rows with the same length as mean and as - # many rows are necessary to form a matrix of shape final_shape. - x = self.standard_normal(np.multiply.reduce(final_shape)) - x.shape = (np.multiply.reduce(final_shape[0:len(final_shape)-1]), - mean.shape[0]) + x = self.standard_normal(final_shape).reshape(-1, mean.shape[0]) + # Transform matrix of standard normals into matrix where each row # contains multivariate normals with the desired covariance. # Compute A such that dot(transpose(A),A) == cov. # Then the matrix products of the rows of x and A has the desired # covariance. Note that sqrt(s)*v where (u,s,v) is the singular value # decomposition of cov is such an A. - - from numpy.dual import svd - # XXX: we really should be doing this by Cholesky decomposition - (u,s,v) = svd(cov) - x = np.dot(x*np.sqrt(s),v) - # The rows of x now have the correct covariance but mean 0. Add - # mean to each row. Then each row will have mean mean. - np.add(mean,x,x) + # + # Also check that cov is positive-semidefinite. If so, the u.T and v + # matrices should be equal up to roundoff error if cov is + # symmetrical and the singular value of the corresponding row is + # not zero. We continue to use the SVD rather than Cholesky in + # order to preserve current outputs. Note that symmetry has not + # been checked. + (u, s, v) = svd(cov) + neg = (np.sum(u.T * v, axis=1) < 0) & (s > 0) + if np.any(neg): + s[neg] = 0. + warnings.warn("covariance is not positive-semidefinite.", + RuntimeWarning) + + x = np.dot(x, np.sqrt(s)[:, None] * v) + x += mean x.shape = tuple(final_shape) return x @@ -4235,7 +4390,7 @@ cdef class RandomState: cdef ndarray parr "arrayObject_parr", mnarr "arrayObject_mnarr" cdef double *pix cdef long *mnix - cdef npy_intp i, j, dn + cdef npy_intp i, j, dn, sz cdef double Sum d = len(pvals) @@ -4245,30 +4400,27 @@ cdef class RandomState: if kahan_sum(pix, d-1) > (1.0 + 1e-12): raise ValueError("sum(pvals[:-1]) > 1.0") - if size is None: - shape = (d,) - elif type(size) is int: - shape = (size, d) - else: - shape = size + (d,) + shape = _shape_from_size(size, d) multin = np.zeros(shape, int) mnarr = <ndarray>multin mnix = <long*>PyArray_DATA(mnarr) - i = 0 - while i < PyArray_SIZE(mnarr): - Sum = 1.0 - dn = n - for j from 0 <= j < d-1: - mnix[i+j] = rk_binomial(self.internal_state, dn, pix[j]/Sum) - dn = dn - mnix[i+j] - if dn <= 0: - break - Sum = Sum - pix[j] - if dn > 0: - mnix[i+d-1] = dn - - i = i + d + sz = PyArray_SIZE(mnarr) + with self.lock, nogil: + i = 0 + while i < sz: + Sum = 1.0 + dn = n + for j from 0 <= j < d-1: + mnix[i+j] = rk_binomial(self.internal_state, dn, pix[j]/Sum) + dn = dn - mnix[i+j] + if dn <= 0: + break + Sum = Sum - pix[j] + if dn > 0: + mnix[i+d-1] = dn + + i = i + d return multin @@ -4354,7 +4506,8 @@ cdef class RandomState: cdef npy_intp k cdef npy_intp totsize cdef ndarray alpha_arr, val_arr - cdef double *alpha_data, *val_data + cdef double *alpha_data + cdef double *val_data cdef npy_intp i, j cdef double acc, invacc @@ -4362,12 +4515,7 @@ cdef class RandomState: alpha_arr = <ndarray>PyArray_ContiguousFromObject(alpha, NPY_DOUBLE, 1, 1) alpha_data = <double*>PyArray_DATA(alpha_arr) - if size is None: - shape = (k,) - elif type(size) is int: - shape = (size, k) - else: - shape = size + (k,) + shape = _shape_from_size(size, k) diric = np.zeros(shape, np.float64) val_arr = <ndarray>diric @@ -4375,15 +4523,17 @@ cdef class RandomState: i = 0 totsize = PyArray_SIZE(val_arr) - while i < totsize: - acc = 0.0 - for j from 0 <= j < k: - val_data[i+j] = rk_standard_gamma(self.internal_state, alpha_data[j]) - acc = acc + val_data[i+j] - invacc = 1/acc - for j from 0 <= j < k: - val_data[i+j] = val_data[i+j] * invacc - i = i + k + with self.lock, nogil: + while i < totsize: + acc = 0.0 + for j from 0 <= j < k: + val_data[i+j] = rk_standard_gamma(self.internal_state, + alpha_data[j]) + acc = acc + val_data[i+j] + invacc = 1/acc + for j from 0 <= j < k: + val_data[i+j] = val_data[i+j] * invacc + i = i + k return diric @@ -4426,11 +4576,12 @@ cdef class RandomState: i = len(x) - 1 # Logic adapted from random.shuffle() - if isinstance(x, np.ndarray) and x.ndim > 1: + if isinstance(x, np.ndarray) and \ + (x.ndim > 1 or x.dtype.fields is not None): # For a multi-dimensional ndarray, indexing returns a view onto # each row. So we can't just use ordinary assignment to swap the # rows; we need a bounce buffer. - buf = np.empty(x.shape[1:], dtype=x.dtype) + buf = np.empty_like(x[0]) while i > 0: j = rk_interval(i, self.internal_state) buf[...] = x[j] |