diff options
Diffstat (limited to 'numpy')
| -rw-r--r-- | numpy/random/include/legacy-distributions.h | 1 | ||||
| -rw-r--r-- | numpy/random/mtrand.pyx | 4 | ||||
| -rw-r--r-- | numpy/random/src/distributions/distributions.c | 29 | ||||
| -rw-r--r-- | numpy/random/src/legacy/legacy-distributions.c | 68 |
4 files changed, 82 insertions, 20 deletions
diff --git a/numpy/random/include/legacy-distributions.h b/numpy/random/include/legacy-distributions.h index b8ba0841c..f7ccd2cb5 100644 --- a/numpy/random/include/legacy-distributions.h +++ b/numpy/random/include/legacy-distributions.h @@ -31,6 +31,7 @@ extern double legacy_f(aug_bitgen_t *aug_state, double dfnum, double dfden); extern double legacy_normal(aug_bitgen_t *aug_state, double loc, double scale); extern double legacy_standard_gamma(aug_bitgen_t *aug_state, double shape); extern double legacy_exponential(aug_bitgen_t *aug_state, double scale); +extern double legacy_vonmises(bitgen_t *bitgen_state, double mu, double kappa); extern int64_t legacy_random_binomial(bitgen_t *bitgen_state, double p, int64_t n, binomial_t *binomial); extern int64_t legacy_negative_binomial(aug_bitgen_t *aug_state, double n, diff --git a/numpy/random/mtrand.pyx b/numpy/random/mtrand.pyx index 8d349b7d8..6f44e271f 100644 --- a/numpy/random/mtrand.pyx +++ b/numpy/random/mtrand.pyx @@ -51,7 +51,6 @@ cdef extern from "numpy/random/distributions.h": void random_standard_uniform_fill(bitgen_t* bitgen_state, np.npy_intp cnt, double *out) nogil int64_t random_positive_int(bitgen_t *bitgen_state) nogil double random_uniform(bitgen_t *bitgen_state, double lower, double range) nogil - double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) nogil double random_laplace(bitgen_t *bitgen_state, double loc, double scale) nogil double random_gumbel(bitgen_t *bitgen_state, double loc, double scale) nogil double random_logistic(bitgen_t *bitgen_state, double loc, double scale) nogil @@ -100,6 +99,7 @@ cdef extern from "include/legacy-distributions.h": double legacy_f(aug_bitgen_t *aug_state, double dfnum, double dfden) nogil double legacy_exponential(aug_bitgen_t *aug_state, double scale) nogil double legacy_power(aug_bitgen_t *state, double a) nogil + double legacy_vonmises(bitgen_t *bitgen_state, double mu, double kappa) nogil np.import_array() @@ -2281,7 +2281,7 @@ cdef class RandomState: >>> plt.show() """ - return cont(&random_vonmises, &self._bitgen, size, self.lock, 2, + return cont(&legacy_vonmises, &self._bitgen, size, self.lock, 2, mu, 'mu', CONS_NONE, kappa, 'kappa', CONS_NON_NEGATIVE, 0.0, '', CONS_NONE, None) diff --git a/numpy/random/src/distributions/distributions.c b/numpy/random/src/distributions/distributions.c index ce909933e..f47c54a53 100644 --- a/numpy/random/src/distributions/distributions.c +++ b/numpy/random/src/distributions/distributions.c @@ -1,7 +1,6 @@ #include "numpy/random/distributions.h" #include "ziggurat_constants.h" #include "logfactorial.h" -#include <stdio.h> #if defined(_MSC_VER) && defined(_WIN64) #include <intrin.h> @@ -844,6 +843,7 @@ double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) { return NPY_NAN; } if (kappa < 1e-8) { + /* Use a uniform for very small values of kappa */ return M_PI * (2 * next_double(bitgen_state) - 1); } else { /* with double precision rho is zero until 1.4e-8 */ @@ -854,30 +854,23 @@ double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) { */ s = (1. / kappa + kappa); } else { -#ifndef NP_RANDOM_LEGACY - /* Fallback to normal distribution for big values of kappa - * Fix only applies to Generator, not RandomState - * RandomState defines NP_RANDOM_LEGACY=1 - */ - if (kappa > 1e6) { + if (kappa <= 1e6) { + /* Path for 1e-5 <= kappa <= 1e6 */ + double r = 1 + sqrt(1 + 4 * kappa * kappa); + double rho = (r - sqrt(2 * r)) / (2 * kappa); + s = (1 + rho * rho) / (2 * rho); + } else { + /* Fallback to wrapped normal distribution for kappa > 1e6 */ result = mu + sqrt(1. / kappa) * random_standard_normal(bitgen_state); - /* Check if result is within bounds */ + /* Ensure result is within bounds */ if (result < -M_PI) { - return result + 2*M_PI; + result += 2*M_PI; } if (result > M_PI) { - return result - 2*M_PI; + result -= 2*M_PI; } return result; - } else { -#endif - /* Path for 1e-5 <= kappa <= 1e6 */ - double r = 1 + sqrt(1 + 4 * kappa * kappa); - double rho = (r - sqrt(2 * r)) / (2 * kappa); - s = (1 + rho * rho) / (2 * rho); -#if !defined NP_RANDOM_LEGACY } -#endif } while (1) { diff --git a/numpy/random/src/legacy/legacy-distributions.c b/numpy/random/src/legacy/legacy-distributions.c index fd067fe8d..5b17401dd 100644 --- a/numpy/random/src/legacy/legacy-distributions.c +++ b/numpy/random/src/legacy/legacy-distributions.c @@ -1,3 +1,13 @@ +/* + * This file contains generation code for distribution that have been modified + * since Generator was introduced. These are preserved using identical code + * to what was in NumPy 1.16 so that the stream of values generated by + * RandomState is not changed when there are changes that affect Generator. + * + * These functions should not be changed except if they contain code that + * cannot be compiled. They should not be changed for bug fixes, performance + * improvements that can change the values produced, or enhancements to precision. + */ #include "include/legacy-distributions.h" @@ -390,3 +400,61 @@ int64_t legacy_random_logseries(bitgen_t *bitgen_state, double p) { binomial_t *binomial) { return random_multinomial(bitgen_state, n, mnix, pix, d, binomial); } + +double legacy_vonmises(bitgen_t *bitgen_state, double mu, double kappa) { + double s; + double U, V, W, Y, Z; + double result, mod; + int neg; + if (npy_isnan(kappa)) { + return NPY_NAN; + } + if (kappa < 1e-8) { + return M_PI * (2 * next_double(bitgen_state) - 1); + } else { + /* with double precision rho is zero until 1.4e-8 */ + if (kappa < 1e-5) { + /* + * second order taylor expansion around kappa = 0 + * precise until relatively large kappas as second order is 0 + */ + s = (1. / kappa + kappa); + } else { + /* Path for 1e-5 <= kappa <= 1e6 */ + double r = 1 + sqrt(1 + 4 * kappa * kappa); + double rho = (r - sqrt(2 * r)) / (2 * kappa); + s = (1 + rho * rho) / (2 * rho); + } + + while (1) { + U = next_double(bitgen_state); + Z = cos(M_PI * U); + W = (1 + s * Z) / (s + Z); + Y = kappa * (s - W); + V = next_double(bitgen_state); + /* + * V==0.0 is ok here since Y >= 0 always leads + * to accept, while Y < 0 always rejects + */ + if ((Y * (2 - Y) - V >= 0) || (log(Y / V) + 1 - Y >= 0)) { + break; + } + } + + U = next_double(bitgen_state); + + result = acos(W); + if (U < 0.5) { + result = -result; + } + result += mu; + neg = (result < 0); + mod = fabs(result); + mod = (fmod(mod + M_PI, 2 * M_PI) - M_PI); + if (neg) { + mod *= -1; + } + + return mod; + } +}
\ No newline at end of file |