From adf717874f6eab7b19a9fa5b2d338dfb63847077 Mon Sep 17 00:00:00 2001 From: Nikita Grigorian Date: Wed, 4 Nov 2026 11:46:55 -0800 Subject: [PATCH] Add the continuous distributions to dpnp.random.Generator Port NumPy's continuous samplers to device kernels: exponential, gamma, beta, chisquare, standard_t, f, lognormal, standard_cauchy, wald, pareto, weibull, power, rayleigh, laplace, gumbel, logistic, triangular and vonmises. The rejection samplers return NaN for device parameters that would make them loop forever. --- CHANGELOG.md | 1 + doc/reference/random.rst | 24 +- dpnp/libtensor/AGENTS.md | 2 +- dpnp/libtensor/extensions/rng/CMakeLists.txt | 2 + .../extensions/rng/distributions.hpp | 2 + dpnp/libtensor/extensions/rng/exponential.cpp | 65 +++ dpnp/libtensor/extensions/rng/gamma.cpp | 69 +++ .../extensions/rng/kernels/distributions.hpp | 297 +++++++++++ .../extensions/rng/kernels/samplers.hpp | 384 ++++++++++++++ dpnp/libtensor/extensions/rng/normal.cpp | 18 + dpnp/libtensor/extensions/rng/rng_py.cpp | 2 + dpnp/libtensor/extensions/rng/uniform.cpp | 35 ++ dpnp/random/_generator.py | 469 ++++++++++++++++++ dpnp/tests/test_random_generator.py | 229 +++++++++ dpnp/tests/test_sycl_queue.py | 22 + dpnp/tests/test_usm_type.py | 2 + 16 files changed, 1621 insertions(+), 2 deletions(-) create mode 100644 dpnp/libtensor/extensions/rng/exponential.cpp create mode 100644 dpnp/libtensor/extensions/rng/gamma.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index e811916e1b4..596f483596a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * Added support for boolean input in `dpnp.positive` and the `+` operator of `dpnp.ndarray` [#3128](https://github.com/IntelPython/dpnp/pull/3128) * Added `dpnp.top_k` to return the `k` largest or smallest elements and their indices along an axis, with a `sorted` keyword to return them in sorted order [#3133](https://github.com/IntelPython/dpnp/pull/3133) * Added the `dpnp.random.Generator` API with `default_rng`, `BitGenerator`, `SeedSequence` and the `Philox4x3210` bit generator, sampling on the device with the oneMKL device RNG API [#3157](https://github.com/IntelPython/dpnp/pull/3157) +* Added the continuous distributions of `numpy.random.Generator`, such as `gamma`, `beta` and `vonmises`, to `dpnp.random.Generator` [#3158](https://github.com/IntelPython/dpnp/pull/3158) ### Changed diff --git a/doc/reference/random.rst b/doc/reference/random.rst index d96b52654d2..69cf865d7fb 100644 --- a/doc/reference/random.rst +++ b/doc/reference/random.rst @@ -37,7 +37,9 @@ Differences from NumPy * Parameters on the host (Python scalars, lists and NumPy arrays) are checked as in NumPy. Parameters in :class:`dpnp.ndarray` are not, since that would wait for the device. - Invalid ones give undefined values, such as NaN. + Invalid ones give undefined values, such as NaN or the smallest ``int64``, + but the kernels handle them safely: they never loop without end or access + memory out of bounds. Random Generator ---------------- @@ -65,9 +67,29 @@ Distributions :toctree: generated/ :nosignatures: + Generator.beta + Generator.chisquare + Generator.exponential + Generator.f + Generator.gamma + Generator.gumbel + Generator.laplace + Generator.logistic + Generator.lognormal Generator.normal + Generator.pareto + Generator.power + Generator.rayleigh + Generator.standard_cauchy + Generator.standard_exponential + Generator.standard_gamma Generator.standard_normal + Generator.standard_t + Generator.triangular Generator.uniform + Generator.vonmises + Generator.wald + Generator.weibull Bit generators -------------- diff --git a/dpnp/libtensor/AGENTS.md b/dpnp/libtensor/AGENTS.md index 556cd6bf62f..0bb184f73f9 100644 --- a/dpnp/libtensor/AGENTS.md +++ b/dpnp/libtensor/AGENTS.md @@ -23,7 +23,7 @@ A new extension's `CMakeLists.txt` lists its sources and calls `build_dpnp_pybin - Guard double, complex and half with `dpnp::tensor::type_utils::validate_type_for_device(q)`; don't use double literals or casts in kernels that must run without fp64. - `blas`, `fft`, `lapack` and `sparse` build against either oneMKL or oneMath (`USE_ONEMATH`); `vm` has kernels only with oneMKL: with oneMath it builds empty and `_is_available()` returns `False`; `rng` uses the header-only oneMKL device RNG API and links no oneMKL library. - Release the GIL around blocking oneMKL calls (`py::gil_scoped_release`). -- `rng` values depend only on the engine state and the number of elements, never on the device or the layout chosen for it (`extensions/rng/kernels/sample.hpp`); each engine's `items_per_work_item` (`kernels/engines.hpp`) is part of its stream definition. +- `rng` values depend only on the engine state and the number of elements, never on the device or the layout chosen for it (`extensions/rng/kernels/sample.hpp`); each engine's `items_per_work_item` (`kernels/engines.hpp`) is part of its stream definition. Parameters on the device are unchecked, so a rejection sampler returns NaN for the inputs that would make it loop forever, including values that are valid but round badly in fp32. - An extension whose kernels may be launched over more than `INT_MAX` work-items passes `LARGE_RANGES` to `build_dpnp_pybind11_ext`, which adds `-fno-sycl-id-queries-fit-in-int` at compile and link time (DPC++ 2026.2 and later need it at link time). `source/CMakeLists.txt` and `dpnp_backend_c` pass it to both `target_compile_options` and `target_link_options`. - clang-format 23 formats C++ (`.clang-format`: LLVM base, 4-space indent, 80 columns); gersemi formats CMake at 88 columns, except `dpnp/backend/cmake/Modules/`. diff --git a/dpnp/libtensor/extensions/rng/CMakeLists.txt b/dpnp/libtensor/extensions/rng/CMakeLists.txt index 20f793c978d..7e41305df9d 100644 --- a/dpnp/libtensor/extensions/rng/CMakeLists.txt +++ b/dpnp/libtensor/extensions/rng/CMakeLists.txt @@ -29,6 +29,8 @@ set(python_module_name _rng_impl) set(_module_src + ${CMAKE_CURRENT_SOURCE_DIR}/exponential.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/gamma.cpp ${CMAKE_CURRENT_SOURCE_DIR}/normal.cpp ${CMAKE_CURRENT_SOURCE_DIR}/rng_py.cpp ${CMAKE_CURRENT_SOURCE_DIR}/sample.cpp diff --git a/dpnp/libtensor/extensions/rng/distributions.hpp b/dpnp/libtensor/extensions/rng/distributions.hpp index 381534f6fe4..6d41e71b39d 100644 --- a/dpnp/libtensor/extensions/rng/distributions.hpp +++ b/dpnp/libtensor/extensions/rng/distributions.hpp @@ -36,4 +36,6 @@ namespace py = pybind11; void init_uniform(py::module_ m); void init_normal(py::module_ m); +void init_exponential(py::module_ m); +void init_gamma(py::module_ m); } // namespace dpnp::extensions::rng diff --git a/dpnp/libtensor/extensions/rng/exponential.cpp b/dpnp/libtensor/extensions/rng/exponential.cpp new file mode 100644 index 00000000000..5e36e15e229 --- /dev/null +++ b/dpnp/libtensor/extensions/rng/exponential.cpp @@ -0,0 +1,65 @@ +//***************************************************************************** +// Copyright (c) 2026, Intel Corporation +// All rights reserved. +// +// 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 the copyright holder 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 HOLDER 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. +//***************************************************************************** +#include + +#include "common.hpp" +#include "distributions.hpp" + +namespace dpnp::extensions::rng +{ +static sample_dispatch_table_t exponential_dispatch_table; +static sample_dispatch_table_t pareto_dispatch_table; +static sample_dispatch_table_t weibull_dispatch_table; +static sample_dispatch_table_t power_dispatch_table; +static sample_dispatch_table_t rayleigh_dispatch_table; + +void init_exponential(py::module_ m) +{ + def_sampler( + m, "_exponential", + "Fill the result with exponential samples of scale scale", + exponential_dispatch_table); + + def_sampler( + m, "_pareto", "Fill the result with Pareto II samples of shape a", + pareto_dispatch_table); + + def_sampler( + m, "_weibull", "Fill the result with Weibull samples of shape a", + weibull_dispatch_table); + + def_sampler( + m, "_power", + "Fill the result with power distribution samples of exponent a - 1", + power_dispatch_table); + + def_sampler( + m, "_rayleigh", "Fill the result with Rayleigh samples of mode scale", + rayleigh_dispatch_table); +} +} // namespace dpnp::extensions::rng diff --git a/dpnp/libtensor/extensions/rng/gamma.cpp b/dpnp/libtensor/extensions/rng/gamma.cpp new file mode 100644 index 00000000000..b683f78bdc2 --- /dev/null +++ b/dpnp/libtensor/extensions/rng/gamma.cpp @@ -0,0 +1,69 @@ +//***************************************************************************** +// Copyright (c) 2026, Intel Corporation +// All rights reserved. +// +// 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 the copyright holder 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 HOLDER 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. +//***************************************************************************** +#include + +#include "common.hpp" +#include "distributions.hpp" + +namespace dpnp::extensions::rng +{ +static sample_dispatch_table_t gamma_dispatch_table; +static sample_dispatch_table_t beta_dispatch_table; +static sample_dispatch_table_t chisquare_dispatch_table; +static sample_dispatch_table_t standard_t_dispatch_table; +static sample_dispatch_table_t f_dispatch_table; + +void init_gamma(py::module_ m) +{ + def_sampler( + m, "_gamma", + "Fill the result with gamma samples of shape shape and scale scale", + gamma_dispatch_table); + + def_sampler( + m, "_beta", + "Fill the result with beta samples of shape parameters a and b", + beta_dispatch_table); + + def_sampler( + m, "_chisquare", + "Fill the result with chi-square samples of df degrees of freedom", + chisquare_dispatch_table); + + def_sampler( + m, "_standard_t", + "Fill the result with Student's t samples of df degrees of freedom", + standard_t_dispatch_table); + + def_sampler( + m, "_f", + "Fill the result with F samples of dfnum and dfden degrees of " + "freedom", + f_dispatch_table); +} +} // namespace dpnp::extensions::rng diff --git a/dpnp/libtensor/extensions/rng/kernels/distributions.hpp b/dpnp/libtensor/extensions/rng/kernels/distributions.hpp index 77bcdd9300b..a39bf73e25e 100644 --- a/dpnp/libtensor/extensions/rng/kernels/distributions.hpp +++ b/dpnp/libtensor/extensions/rng/kernels/distributions.hpp @@ -64,6 +64,7 @@ #pragma once #include +#include #include @@ -111,4 +112,300 @@ class StandardNormal return r * sycl::sinpi(phi); } }; + +// NumPy's "inv" method; 1 - u is on (0, 1], so the result is finite +template +T standard_exponential(BitGenT &bg) +{ + return -sycl::log1p(-bg.template next_uniform()); +} + +// Marsaglia and Tsang for shape > 1, and a rejection method built on the +// exponential distribution for shape < 1, as in NumPy's +// random_standard_gamma() +template +T standard_gamma(BitGenT &bg, StandardNormal &normal, const T shape) +{ + if (shape == T(1)) { + return standard_exponential(bg); + } + else if (shape == T(0)) { + return T(0); + } + else if (!(shape > T(0))) { + // invalid, and for a shape just below zero every draw would be + // rejected + return std::numeric_limits::quiet_NaN(); + } + else if (shape < T(1)) { + for (;;) { + const T U = bg.template next_uniform(); + const T V = standard_exponential(bg); + if (U <= T(1) - shape) { + const T X = sycl::pow(U, T(1) / shape); + if (X <= V) { + return X; + } + } + else { + const T Y = -sycl::log((T(1) - U) / shape); + const T X = sycl::pow(T(1) - shape + shape * Y, T(1) / shape); + if (X <= (V + Y)) { + return X; + } + } + } + } + + const T b = shape - T(1) / T(3); + const T c = T(1) / sycl::sqrt(T(9) * b); + for (;;) { + T X, V; + do { + X = normal(bg); + V = T(1) + c * X; + } while (V <= T(0)); + + V = V * V * V; + const T U = bg.template next_uniform(); + if (U < T(1) - T(0.0331) * (X * X) * (X * X)) { + return b * V; + } + // log(0) is fine here + if (sycl::log(U) < T(0.5) * X * X + b * (T(1) - V + sycl::log(V))) { + return b * V; + } + } +} + +// Johnk's algorithm when both parameters are at most one, and a ratio of +// gamma variates otherwise, as in NumPy's random_beta() +template +T beta(BitGenT &bg, StandardNormal &normal, const T a, const T b) +{ + if (!(a > T(0)) || !(b > T(0))) { + // invalid input, which Johnk's algorithm would reject forever + return std::numeric_limits::quiet_NaN(); + } + if ((a <= T(1)) && (b <= T(1))) { + // Below this threshold the result is 0 or 1 with overwhelming + // probability, so a single uniform draw decides which. The + // threshold underflows to zero, disabling the shortcut, in float. + constexpr T tiny_threshold = static_cast(3e-103); + if (a < tiny_threshold && b < tiny_threshold) { + const T U = bg.template next_uniform(); + return static_cast((a + b) * U < a); + } + + for (;;) { + const T U = bg.template next_uniform(); + const T V = bg.template next_uniform(); + const T X = sycl::pow(U, T(1) / a); + const T Y = sycl::pow(V, T(1) / b); + const T XpY = X + Y; + // reject if both U and V are 0 + if ((XpY <= T(1)) && (U + V > T(0))) { + T res; + + // NumPy tests X > 0 and Y > 0 here. Devices that flush + // subnormals to zero can see a subnormal X as positive + // while X + Y flushes to zero, so require normal values. + constexpr T min_normal = std::numeric_limits::min(); + if ((X >= min_normal) && (Y >= min_normal)) { + res = X / XpY; + } + else { + // X or Y underflowed, so work with logarithms + const T logX = sycl::log(U) / a; + const T logY = sycl::log(V) / b; + const T delta = logX - logY; + if (delta > T(0)) { + res = sycl::exp(-sycl::log1p(sycl::exp(-delta))); + } + else { + res = sycl::exp(delta - sycl::log1p(sycl::exp(delta))); + } + } + + // Division and exp need not be correctly rounded on the + // device, and can land an ulp above 1; NaN is kept + return (res > T(1)) ? T(1) : res; + } + } + } + + const T Ga = standard_gamma(bg, normal, a); + const T Gb = standard_gamma(bg, normal, b); + // as above; this is Ga / Ga when Gb underflows + const T res = Ga / (Ga + Gb); + return (res > T(1)) ? T(1) : res; +} + +template +T chisquare(BitGenT &bg, StandardNormal &normal, const T df) +{ + return T(2) * standard_gamma(bg, normal, df / T(2)); +} + +template +T standard_cauchy(BitGenT &bg, StandardNormal &normal) +{ + const T num = normal(bg); + return num / normal(bg); +} + +template +T standard_t(BitGenT &bg, StandardNormal &normal, const T df) +{ + const T num = normal(bg); + const T denom = standard_gamma(bg, normal, df / T(2)); + return sycl::sqrt(df / T(2)) * num / sycl::sqrt(denom); +} + +template +T f(BitGenT &bg, StandardNormal &normal, const T dfnum, const T dfden) +{ + const T num = chisquare(bg, normal, dfnum) * dfden; + return num / (chisquare(bg, normal, dfden) * dfnum); +} + +template +T laplace(BitGenT &bg, const T loc, const T scale) +{ + for (;;) { + const T U = bg.template next_uniform(); + if (U >= T(0.5)) { + return loc - scale * sycl::log(T(2) - U - U); + } + else if (U > T(0)) { + return loc + scale * sycl::log(U + U); + } + // reject U == 0 + } +} + +template +T gumbel(BitGenT &bg, const T loc, const T scale) +{ + for (;;) { + const T U = T(1) - bg.template next_uniform(); + if (U < T(1)) { + return loc - scale * sycl::log(-sycl::log(U)); + } + // reject U == 1 + } +} + +template +T logistic(BitGenT &bg, const T loc, const T scale) +{ + for (;;) { + const T U = bg.template next_uniform(); + if (U > T(0)) { + return loc + scale * sycl::log(U / (T(1) - U)); + } + // reject U == 0 + } +} + +// NumPy's random_wald() takes X = mean + mean/(2*scale)*(Y - sqrt(4*scale*Y + +// Y*Y)), which cancels, and overflows in float, for large Y. This is the same +// value as mean*4*scale/(sqrt(Y) + sqrt(Y + 4*scale))**2. +template +T wald(BitGenT &bg, StandardNormal &normal, const T mean, const T scale) +{ + const T n = normal(bg); + const T Y = mean * n * n; + const T s = sycl::sqrt(Y) + sycl::sqrt(Y + T(4) * scale); + const T X = mean * (T(4) * scale / (s * s)); + const T U = bg.template next_uniform(); + if (U <= mean / (mean + X)) { + return X; + } + return mean * (mean / X); +} + +template +T triangular(BitGenT &bg, const T left, const T mode, const T right) +{ + const T base = right - left; + const T leftbase = mode - left; + const T ratio = leftbase / base; + const T leftprod = leftbase * base; + const T rightprod = (right - mode) * base; + + const T U = bg.template next_uniform(); + if (U <= ratio) { + return left + sycl::sqrt(U * leftprod); + } + return right - sycl::sqrt((T(1) - U) * rightprod); +} + +// Best and Fisher's algorithm, with the fallbacks of NumPy's +// random_vonmises() for very small and very large kappa +template +T vonmises(BitGenT &bg, StandardNormal &normal, const T mu, const T kappa) +{ + constexpr T pi = T(3.141592653589793238462643383279502884L); + + if (sycl::isnan(kappa)) { + return std::numeric_limits::quiet_NaN(); + } + if (kappa < T(1e-8)) { + // the distribution is uniform on the circle + return pi * (T(2) * bg.template next_uniform() - T(1)); + } + + T s; + if (kappa < T(1e-5)) { + // as in NumPy, whose form of rho is zero in double precision until + // 1.4e-8 + s = T(1) / kappa + kappa; + } + else if (kappa <= T(1e6)) { + // rho = (r - sqrt(2 r)) / (2 kappa), without the cancellation which + // makes it zero in single precision for kappa below about 4e-4, + // where every draw would be rejected + const T r = T(1) + sycl::sqrt(T(1) + T(4) * kappa * kappa); + const T rho = T(2) * kappa / (r + sycl::sqrt(T(2) * r)); + s = (T(1) + rho * rho) / (T(2) * rho); + } + else { + // a wrapped normal approximates the distribution for large kappa + T result = mu + sycl::sqrt(T(1) / kappa) * normal(bg); + if (result < -pi) { + result += T(2) * pi; + } + if (result > pi) { + result -= T(2) * pi; + } + return result; + } + + T W; + for (;;) { + const T U = bg.template next_uniform(); + const T Z = sycl::cospi(U); + W = (T(1) + s * Z) / (s + Z); + const T Y = kappa * (s - W); + const T V = bg.template next_uniform(); + // V == 0 is fine: Y >= 0 always accepts and Y < 0 always rejects + if ((Y * (T(2) - Y) - V >= T(0)) || + (sycl::log(Y / V) + T(1) - Y >= T(0))) { + break; + } + } + + const T U = bg.template next_uniform(); + T result = sycl::acos(W); + if (U < T(0.5)) { + result = -result; + } + result += mu; + + const bool neg = (result < T(0)); + T mod = sycl::fabs(result); + mod = sycl::fmod(mod + pi, T(2) * pi) - pi; + return neg ? -mod : mod; +} } // namespace dpnp::kernels::rng diff --git a/dpnp/libtensor/extensions/rng/kernels/samplers.hpp b/dpnp/libtensor/extensions/rng/kernels/samplers.hpp index bd66f2bb3ae..c63cf4f2b8d 100644 --- a/dpnp/libtensor/extensions/rng/kernels/samplers.hpp +++ b/dpnp/libtensor/extensions/rng/kernels/samplers.hpp @@ -131,4 +131,388 @@ class NormalSampler : public SamplerBase return loc_[i] + scale_[i] * normal_(bg); } }; + +// Log-normal: exp of a normal with mean `mean` and standard deviation `sigma` +template +class LognormalSampler : public SamplerBase +{ +private: + Param mean_; + Param sigma_; + StandardNormal normal_{}; + +public: + // one normal, half a pair of uniform values, per element + static constexpr unsigned words_per_item = uniform_words; + + LognormalSampler(const Param &mean, const Param &sigma) + : mean_(mean), sigma_(sigma) + { + } + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return sycl::exp(mean_[i] + sigma_[i] * normal_(bg)); + } +}; + +template +class StandardCauchySampler : public SamplerBase +{ +private: + StandardNormal normal_{}; + +public: + // two normals per element + static constexpr unsigned words_per_item = 2 * uniform_words; + + template + T operator()(BitGenT &bg, const std::size_t) + { + return standard_cauchy(bg, normal_); + } +}; + +// Wald (inverse Gaussian) with mean `mean` and shape `scale` +template +class WaldSampler : public SamplerBase +{ +private: + Param mean_; + Param scale_; + StandardNormal normal_{}; + +public: + // one normal and one uniform value per element + static constexpr unsigned words_per_item = 2 * uniform_words; + + WaldSampler(const Param &mean, const Param &scale) + : mean_(mean), scale_(scale) + { + } + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return wald(bg, normal_, mean_[i], scale_[i]); + } +}; + +// Exponential with scale `scale` +template +class ExponentialSampler : public SamplerBase +{ +private: + Param scale_; + +public: + // one uniform value per element + static constexpr unsigned words_per_item = uniform_words; + + explicit ExponentialSampler(const Param &scale) : scale_(scale) {} + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return scale_[i] * standard_exponential(bg); + } +}; + +// Pareto II (Lomax) with shape `a` +template +class ParetoSampler : public SamplerBase +{ +private: + Param a_; + +public: + // one uniform value per element + static constexpr unsigned words_per_item = uniform_words; + + explicit ParetoSampler(const Param &a) : a_(a) {} + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return sycl::expm1(standard_exponential(bg) / a_[i]); + } +}; + +// Weibull with shape `a` +template +class WeibullSampler : public SamplerBase +{ +private: + Param a_; + +public: + // one uniform value per element + static constexpr unsigned words_per_item = uniform_words; + + explicit WeibullSampler(const Param &a) : a_(a) {} + + template + T operator()(BitGenT &bg, const std::size_t i) + { + const T a = a_[i]; + const T E = standard_exponential(bg); + if (a == T(0)) { + return T(0); + } + return sycl::pow(E, T(1) / a); + } +}; + +// Power distribution with exponent `a` - 1, on [0, 1] +template +class PowerSampler : public SamplerBase +{ +private: + Param a_; + +public: + // one uniform value per element + static constexpr unsigned words_per_item = uniform_words; + + explicit PowerSampler(const Param &a) : a_(a) {} + + template + T operator()(BitGenT &bg, const std::size_t i) + { + const T E = standard_exponential(bg); + return sycl::pow(-sycl::expm1(-E), T(1) / a_[i]); + } +}; + +// Rayleigh with mode `scale` +template +class RayleighSampler : public SamplerBase +{ +private: + Param scale_; + +public: + // one uniform value per element + static constexpr unsigned words_per_item = uniform_words; + + explicit RayleighSampler(const Param &scale) : scale_(scale) {} + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return scale_[i] * sycl::sqrt(T(2) * standard_exponential(bg)); + } +}; + +// Gamma with shape `shape` and scale `scale` +template +class GammaSampler : public SamplerBase +{ +private: + Param shape_; + Param scale_; + StandardNormal normal_{}; + +public: + GammaSampler(const Param &shape, const Param &scale) + : shape_(shape), scale_(scale) + { + } + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return scale_[i] * standard_gamma(bg, normal_, shape_[i]); + } +}; + +// Beta with shape parameters `a` and `b` +template +class BetaSampler : public SamplerBase +{ +private: + Param a_; + Param b_; + StandardNormal normal_{}; + +public: + BetaSampler(const Param &a, const Param &b) : a_(a), b_(b) {} + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return beta(bg, normal_, a_[i], b_[i]); + } +}; + +// Chi-square with `df` degrees of freedom +template +class ChisquareSampler : public SamplerBase +{ +private: + Param df_; + StandardNormal normal_{}; + +public: + explicit ChisquareSampler(const Param &df) : df_(df) {} + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return chisquare(bg, normal_, df_[i]); + } +}; + +// Student's t with `df` degrees of freedom +template +class StandardTSampler : public SamplerBase +{ +private: + Param df_; + StandardNormal normal_{}; + +public: + explicit StandardTSampler(const Param &df) : df_(df) {} + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return standard_t(bg, normal_, df_[i]); + } +}; + +// F with `dfnum` and `dfden` degrees of freedom +template +class FSampler : public SamplerBase +{ +private: + Param dfnum_; + Param dfden_; + StandardNormal normal_{}; + +public: + FSampler(const Param &dfnum, const Param &dfden) + : dfnum_(dfnum), dfden_(dfden) + { + } + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return f(bg, normal_, dfnum_[i], dfden_[i]); + } +}; + +// Laplace with location `loc` and scale `scale` +template +class LaplaceSampler : public SamplerBase +{ +private: + Param loc_; + Param scale_; + +public: + LaplaceSampler(const Param &loc, const Param &scale) + : loc_(loc), scale_(scale) + { + } + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return laplace(bg, loc_[i], scale_[i]); + } +}; + +// Gumbel with location `loc` and scale `scale` +template +class GumbelSampler : public SamplerBase +{ +private: + Param loc_; + Param scale_; + +public: + GumbelSampler(const Param &loc, const Param &scale) + : loc_(loc), scale_(scale) + { + } + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return gumbel(bg, loc_[i], scale_[i]); + } +}; + +// Logistic with location `loc` and scale `scale` +template +class LogisticSampler : public SamplerBase +{ +private: + Param loc_; + Param scale_; + +public: + LogisticSampler(const Param &loc, const Param &scale) + : loc_(loc), scale_(scale) + { + } + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return logistic(bg, loc_[i], scale_[i]); + } +}; + +// Triangular on [left, right] with mode `mode` +template +class TriangularSampler : public SamplerBase +{ +private: + Param left_; + Param mode_; + Param right_; + +public: + // one uniform value per element + static constexpr unsigned words_per_item = uniform_words; + + TriangularSampler(const Param &left, + const Param &mode, + const Param &right) + : left_(left), mode_(mode), right_(right) + { + } + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return triangular(bg, left_[i], mode_[i], right_[i]); + } +}; + +// von Mises with mode `mu` and concentration `kappa`, on [-pi, pi] +template +class VonmisesSampler : public SamplerBase +{ +private: + Param mu_; + Param kappa_; + StandardNormal normal_{}; + +public: + VonmisesSampler(const Param &mu, const Param &kappa) + : mu_(mu), kappa_(kappa) + { + } + + template + T operator()(BitGenT &bg, const std::size_t i) + { + return vonmises(bg, normal_, mu_[i], kappa_[i]); + } +}; } // namespace dpnp::kernels::rng diff --git a/dpnp/libtensor/extensions/rng/normal.cpp b/dpnp/libtensor/extensions/rng/normal.cpp index b2660ae48ad..df49ab1a032 100644 --- a/dpnp/libtensor/extensions/rng/normal.cpp +++ b/dpnp/libtensor/extensions/rng/normal.cpp @@ -33,6 +33,9 @@ namespace dpnp::extensions::rng { static sample_dispatch_table_t normal_dispatch_table; +static sample_dispatch_table_t lognormal_dispatch_table; +static sample_dispatch_table_t standard_cauchy_dispatch_table; +static sample_dispatch_table_t wald_dispatch_table; void init_normal(py::module_ m) { @@ -41,5 +44,20 @@ void init_normal(py::module_ m) "Fill the result with normal samples of mean loc and standard " "deviation scale", normal_dispatch_table); + + def_sampler( + m, "_lognormal", + "Fill the result with log-normal samples whose logarithm has mean " + "mean and standard deviation sigma", + lognormal_dispatch_table); + + def_sampler( + m, "_standard_cauchy", "Fill the result with standard Cauchy samples", + standard_cauchy_dispatch_table); + + def_sampler( + m, "_wald", + "Fill the result with Wald samples of mean mean and shape scale", + wald_dispatch_table); } } // namespace dpnp::extensions::rng diff --git a/dpnp/libtensor/extensions/rng/rng_py.cpp b/dpnp/libtensor/extensions/rng/rng_py.cpp index eba9388c12a..e7f058d9770 100644 --- a/dpnp/libtensor/extensions/rng/rng_py.cpp +++ b/dpnp/libtensor/extensions/rng/rng_py.cpp @@ -56,4 +56,6 @@ PYBIND11_MODULE(_rng_impl, m, py::mod_gil_not_used()) rng_ns::init_uniform(m); rng_ns::init_normal(m); + rng_ns::init_exponential(m); + rng_ns::init_gamma(m); } diff --git a/dpnp/libtensor/extensions/rng/uniform.cpp b/dpnp/libtensor/extensions/rng/uniform.cpp index 709517fe9bb..629edc5c3b6 100644 --- a/dpnp/libtensor/extensions/rng/uniform.cpp +++ b/dpnp/libtensor/extensions/rng/uniform.cpp @@ -34,6 +34,11 @@ namespace dpnp::extensions::rng { static sample_dispatch_table_t uniform_dispatch_table; static sample_dispatch_table_t random_raw_dispatch_table; +static sample_dispatch_table_t laplace_dispatch_table; +static sample_dispatch_table_t gumbel_dispatch_table; +static sample_dispatch_table_t logistic_dispatch_table; +static sample_dispatch_table_t triangular_dispatch_table; +static sample_dispatch_table_t vonmises_dispatch_table; void init_uniform(py::module_ m) { @@ -46,5 +51,35 @@ void init_uniform(py::module_ m) "Fill the result with the raw 64-bit output of the engine, or with " "its 32-bit halves, low half first", random_raw_dispatch_table); + + def_sampler( + m, "_laplace", + "Fill the result with Laplace samples of location loc and scale " + "scale", + laplace_dispatch_table); + + def_sampler( + m, "_gumbel", + "Fill the result with Gumbel samples of location loc and scale " + "scale", + gumbel_dispatch_table); + + def_sampler( + m, "_logistic", + "Fill the result with logistic samples of location loc and scale " + "scale", + logistic_dispatch_table); + + def_sampler( + m, "_triangular", + "Fill the result with triangular samples on [left, right] with mode " + "mode", + triangular_dispatch_table); + + def_sampler( + m, "_vonmises", + "Fill the result with von Mises samples of mode mu and " + "concentration kappa", + vonmises_dispatch_table); } } // namespace dpnp::extensions::rng diff --git a/dpnp/random/_generator.py b/dpnp/random/_generator.py index a8b087a1248..0367d239ace 100644 --- a/dpnp/random/_generator.py +++ b/dpnp/random/_generator.py @@ -37,6 +37,7 @@ # pylint: disable=no-name-in-module # pylint: disable=protected-access +# pylint: disable=too-many-lines import operator @@ -57,6 +58,7 @@ # numpy/random/_common.pyx _CONS_NONE = 0 _CONS_NON_NEGATIVE = 1 +_CONS_POSITIVE = 2 def _is_array(x): @@ -173,8 +175,12 @@ def _check_constraint(value, name, cons, scalar): if cons == _CONS_NON_NEGATIVE: if numpy.any(~numpy.isnan(val) & numpy.signbit(val)): raise ValueError(f"{name} < 0") + elif cons == _CONS_POSITIVE: + if numpy.any(val <= 0): + raise ValueError(f"{name} <= 0") +# pylint: disable=too-many-public-methods class Generator: """ A container for the bit generator and the distribution samplers. @@ -520,6 +526,469 @@ def normal( (scale, "scale", _CONS_NON_NEGATIVE), ) + def lognormal( + self, mean=0.0, sigma=1.0, size=None, dtype=None, *, usm_type="device" + ): + """ + Draw samples from a log-normal distribution. + + For full documentation refer to + :obj:`numpy.random.Generator.lognormal`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._lognormal, + "lognormal", + size, + dtype, + usm_type, + (mean, "mean", _CONS_NONE), + (sigma, "sigma", _CONS_NON_NEGATIVE), + ) + + def standard_cauchy(self, size=None, dtype=None, *, usm_type="device"): + """ + Draw samples from a standard Cauchy distribution with mode = 0. + + For full documentation refer to + :obj:`numpy.random.Generator.standard_cauchy`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._float_sample( + ri._standard_cauchy, (), size, dtype, usm_type, "standard_cauchy" + ) + + def wald(self, mean, scale, size=None, dtype=None, *, usm_type="device"): + """ + Draw samples from a Wald, or inverse Gaussian, distribution. + + For full documentation refer to :obj:`numpy.random.Generator.wald`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._wald, + "wald", + size, + dtype, + usm_type, + (mean, "mean", _CONS_POSITIVE), + (scale, "scale", _CONS_POSITIVE), + ) + + def standard_exponential( + self, + size=None, + dtype=None, + method="zig", + out=None, + *, + usm_type="device", + ): + """ + Draw samples from the standard exponential distribution. + + For full documentation refer to + :obj:`numpy.random.Generator.standard_exponential`. + + Limitations + ----------- + Both methods use inversion; ``"zig"`` is accepted for compatibility. + + """ + + if method not in ("zig", "inv"): + raise ValueError("method must be either 'inv' or 'zig'") + return self._float_sample( + ri._exponential, + (), + size, + dtype, + usm_type, + "standard_exponential", + out=out, + kernel_params=(1.0,), + ) + + def exponential( + self, scale=1.0, size=None, dtype=None, *, usm_type="device" + ): + """ + Draw samples from an exponential distribution. + + For full documentation refer to + :obj:`numpy.random.Generator.exponential`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._exponential, + "exponential", + size, + dtype, + usm_type, + (scale, "scale", _CONS_NON_NEGATIVE), + ) + + def pareto(self, a, size=None, dtype=None, *, usm_type="device"): + """ + Draw samples from a Pareto II (AKA Lomax) distribution with the + specified shape. + + For full documentation refer to :obj:`numpy.random.Generator.pareto`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._pareto, + "pareto", + size, + dtype, + usm_type, + (a, "a", _CONS_POSITIVE), + ) + + def weibull(self, a, size=None, dtype=None, *, usm_type="device"): + """ + Draw samples from a Weibull distribution. + + For full documentation refer to :obj:`numpy.random.Generator.weibull`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._weibull, + "weibull", + size, + dtype, + usm_type, + (a, "a", _CONS_NON_NEGATIVE), + ) + + def power(self, a, size=None, dtype=None, *, usm_type="device"): + """ + Draw samples in [0, 1] from a power distribution with positive + exponent ``a - 1``. + + For full documentation refer to :obj:`numpy.random.Generator.power`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._power, "power", size, dtype, usm_type, (a, "a", _CONS_POSITIVE) + ) + + def rayleigh(self, scale=1.0, size=None, dtype=None, *, usm_type="device"): + """ + Draw samples from a Rayleigh distribution. + + For full documentation refer to + :obj:`numpy.random.Generator.rayleigh`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._rayleigh, + "rayleigh", + size, + dtype, + usm_type, + (scale, "scale", _CONS_NON_NEGATIVE), + ) + + def laplace( + self, loc=0.0, scale=1.0, size=None, dtype=None, *, usm_type="device" + ): + """ + Draw samples from the Laplace or double exponential distribution + with specified location (or mean) and scale (decay). + + For full documentation refer to :obj:`numpy.random.Generator.laplace`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._laplace, + "laplace", + size, + dtype, + usm_type, + (loc, "loc", _CONS_NONE), + (scale, "scale", _CONS_NON_NEGATIVE), + ) + + def gumbel( + self, loc=0.0, scale=1.0, size=None, dtype=None, *, usm_type="device" + ): + """ + Draw samples from a Gumbel distribution. + + For full documentation refer to :obj:`numpy.random.Generator.gumbel`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._gumbel, + "gumbel", + size, + dtype, + usm_type, + (loc, "loc", _CONS_NONE), + (scale, "scale", _CONS_NON_NEGATIVE), + ) + + def logistic( + self, loc=0.0, scale=1.0, size=None, dtype=None, *, usm_type="device" + ): + """ + Draw samples from a logistic distribution. + + For full documentation refer to + :obj:`numpy.random.Generator.logistic`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._logistic, + "logistic", + size, + dtype, + usm_type, + (loc, "loc", _CONS_NONE), + (scale, "scale", _CONS_NON_NEGATIVE), + ) + + def triangular( + self, left, mode, right, size=None, dtype=None, *, usm_type="device" + ): + """ + Draw samples from the triangular distribution over the interval + [left, right]. + + For full documentation refer to + :obj:`numpy.random.Generator.triangular`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + + params = (left, mode, right) + self._check_queues(*params) + xp, (fleft, fmode, fright) = _check_arrays(params, numpy.float64) + if xp is not None: + if xp is numpy and _is_scalar_path(*params): + for p in params: + _check_real(p) + if xp.any(fleft > fmode): + raise ValueError("left > mode") + if xp.any(fmode > fright): + raise ValueError("mode > right") + if xp.any(fleft == fright): + raise ValueError("left == right") + return self._float_sample( + ri._triangular, params, size, dtype, usm_type, "triangular" + ) + + def vonmises(self, mu, kappa, size=None, dtype=None, *, usm_type="device"): + """ + Draw samples from a von Mises distribution. + + For full documentation refer to + :obj:`numpy.random.Generator.vonmises`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._vonmises, + "vonmises", + size, + dtype, + usm_type, + (mu, "mu", _CONS_NONE), + (kappa, "kappa", _CONS_NON_NEGATIVE), + ) + + def standard_gamma( + self, shape, size=None, dtype=None, out=None, *, usm_type="device" + ): + """ + Draw samples from a standard gamma distribution. + + For full documentation refer to + :obj:`numpy.random.Generator.standard_gamma`. + + Limitations + ----------- + `usm_type` is not in NumPy's signature. + + """ + + self._check_queues(shape, out) + _check_constraint( + shape, + "shape", + _CONS_NON_NEGATIVE, + _is_scalar_path(shape), + ) + return self._float_sample( + ri._gamma, + (shape,), + size, + dtype, + usm_type, + "standard_gamma", + out=out, + kernel_params=(shape, 1.0), + ) + + def gamma( + self, shape, scale=1.0, size=None, dtype=None, *, usm_type="device" + ): + """ + Draw samples from a gamma distribution. + + For full documentation refer to :obj:`numpy.random.Generator.gamma`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._gamma, + "gamma", + size, + dtype, + usm_type, + (shape, "shape", _CONS_NON_NEGATIVE), + (scale, "scale", _CONS_NON_NEGATIVE), + ) + + def chisquare(self, df, size=None, dtype=None, *, usm_type="device"): + """ + Draw samples from a chi-square distribution. + + For full documentation refer to + :obj:`numpy.random.Generator.chisquare`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._chisquare, + "chisquare", + size, + dtype, + usm_type, + (df, "df", _CONS_POSITIVE), + ) + + def standard_t(self, df, size=None, dtype=None, *, usm_type="device"): + """ + Draw samples from a standard Student's t distribution with `df` + degrees of freedom. + + For full documentation refer to + :obj:`numpy.random.Generator.standard_t`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._standard_t, + "standard_t", + size, + dtype, + usm_type, + (df, "df", _CONS_POSITIVE), + ) + + def f(self, dfnum, dfden, size=None, dtype=None, *, usm_type="device"): + """ + Draw samples from an F distribution. + + For full documentation refer to :obj:`numpy.random.Generator.f`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._f, + "f", + size, + dtype, + usm_type, + (dfnum, "dfnum", _CONS_POSITIVE), + (dfden, "dfden", _CONS_POSITIVE), + ) + + def beta(self, a, b, size=None, dtype=None, *, usm_type="device"): + """ + Draw samples from a beta distribution. + + For full documentation refer to :obj:`numpy.random.Generator.beta`. + + Limitations + ----------- + `dtype` and `usm_type` are not in NumPy's signature. + + """ + return self._cont( + ri._beta, + "beta", + size, + dtype, + usm_type, + (a, "a", _CONS_POSITIVE), + (b, "b", _CONS_POSITIVE), + ) + def default_rng(seed=None, *, device=None, sycl_queue=None): """ diff --git a/dpnp/tests/test_random_generator.py b/dpnp/tests/test_random_generator.py index e8e83295f0e..fe68a855ec5 100644 --- a/dpnp/tests/test_random_generator.py +++ b/dpnp/tests/test_random_generator.py @@ -1,3 +1,5 @@ +import math + import dpctl import numpy import pytest @@ -183,6 +185,26 @@ class TestFloatMethods: ("uniform", (-1.0, 2.0)), ("standard_normal", ()), ("normal", (1.0, 2.0)), + ("lognormal", (0.0, 0.5)), + ("standard_exponential", ()), + ("exponential", (2.0,)), + ("standard_gamma", (0.5,)), + ("gamma", (3.0, 2.0)), + ("chisquare", (3.0,)), + ("beta", (0.5, 2.0)), + ("standard_cauchy", ()), + ("wald", (1.0, 2.0)), + ("pareto", (3.0,)), + ("weibull", (2.0,)), + ("power", (3.0,)), + ("rayleigh", (2.0,)), + ("laplace", (1.0, 2.0)), + ("gumbel", (1.0, 2.0)), + ("logistic", (1.0, 2.0)), + ("triangular", (0.0, 1.0, 3.0)), + ("vonmises", (0.5, 2.0)), + ("standard_t", (5.0,)), + ("f", (5.0, 7.0)), ] @pytest.mark.parametrize("dtype", float_dtypes) @@ -234,6 +256,47 @@ def test_random_range(self, dtype): ("uniform", (-1.0, 3.0), 1.0, 16 / 12), ("standard_normal", (), 0.0, 1.0), ("normal", (1.0, 2.0), 1.0, 4.0), + ( + "lognormal", + (0.0, 0.5), + numpy.exp(0.125), + (numpy.exp(0.25) - 1) * numpy.exp(0.25), + ), + ("standard_exponential", (), 1.0, 1.0), + ("exponential", (2.0,), 2.0, 4.0), + ("standard_gamma", (0.5,), 0.5, 0.5), + ("standard_gamma", (1.0,), 1.0, 1.0), + ("gamma", (3.0, 2.0), 6.0, 12.0), + ("chisquare", (3.0,), 3.0, 6.0), + ("beta", (0.5, 2.0), 0.2, 8.0 / 175), + ("beta", (0.5, 0.5), 0.5, 0.125), + ("beta", (2.0, 3.0), 0.4, 0.04), + ("wald", (1.0, 2.0), 1.0, 0.5), + ("wald", (3.0, 50.0), 3.0, 27 / 50), + ("pareto", (20.0,), 1 / 19, 20 / (19**2 * 18)), + ( + "weibull", + (2.0,), + math.gamma(1.5), + math.gamma(2.0) - math.gamma(1.5) ** 2, + ), + ("weibull", (0.0,), 0.0, 0.0), + ("power", (3.0,), 0.75, 0.0375), + ("rayleigh", (2.0,), 2 * math.sqrt(math.pi / 2), 2 * (4 - math.pi)), + ("laplace", (1.0, 2.0), 1.0, 8.0), + ( + "gumbel", + (1.0, 2.0), + 1 + 2 * 0.5772156649015329, + math.pi**2 * 4 / 6, + ), + ("logistic", (1.0, 2.0), 1.0, math.pi**2 * 4 / 3), + ("triangular", (0.0, 1.0, 3.0), 4 / 3, 7 / 18), + ("triangular", (0.0, 0.0, 1.0), 1 / 3, 1 / 18), + # uniform on the circle, whatever mu is, as in NumPy + ("vonmises", (0.5, 0.0), 0.0, math.pi**2 / 3), + ("standard_t", (10.0,), 0.0, 1.25), + ("f", (10.0, 30.0), 30 / 28, 2 * 900 * 38 / (10 * 28**2 * 26)), ], ) def test_moments(self, name, params, mean, var, dtype): @@ -246,6 +309,86 @@ def test_moments(self, name, params, mean, var, dtype): assert abs(res.mean() - mean) <= 5 * numpy.sqrt(var / n) assert abs(res.var() - var) <= 5 * var * numpy.sqrt(10 / n) + @pytest.mark.parametrize("dtype", float_dtypes) + def test_standard_cauchy(self, dtype): + res = default_rng(2).standard_cauchy(size=10**6, dtype=dtype) + q1, q2, q3 = numpy.quantile(dpnp.asnumpy(res), [0.25, 0.5, 0.75]) + assert abs(q1 + 1) < 0.01 and abs(q2) < 0.01 and abs(q3 - 1) < 0.01 + + @pytest.mark.parametrize("dtype", float_dtypes) + @pytest.mark.parametrize("kappa", [1e-9, 1e-6, 0.5, 2.0, 50.0, 1e7]) + def test_vonmises(self, kappa, dtype): + n = 10**6 + mu = 3.0 + res = default_rng(2).vonmises(mu, kappa, size=n, dtype=dtype) + res = dpnp.asnumpy(res).astype(numpy.float64) + assert ((res >= -numpy.pi) & (res <= numpy.pi)).all() + + # E[cos(x - mu)] is I1(kappa) / I0(kappa) + t = numpy.linspace(0, numpy.pi, 20001) + w = numpy.exp(kappa * (numpy.cos(t) - 1)) + expected = numpy.trapezoid(numpy.cos(t) * w, t) / numpy.trapezoid(w, t) + assert abs(numpy.cos(res - mu).mean() - expected) < 5 / numpy.sqrt(n) + # the angles are symmetric around mu + assert abs(numpy.sin(res - mu).mean()) < 5 / numpy.sqrt(n) + + @pytest.mark.parametrize("dtype", float_dtypes) + @pytest.mark.parametrize( + "name, params, lo, hi", + [ + ("pareto", (1e-3,), 0, numpy.inf), + ("weibull", (1e-3,), 0, numpy.inf), + ("power", (1e-3,), 0, 1), + ("power", (1e3,), 0, 1), + ("wald", (1e10, 1e-10), 0, numpy.inf), + ("wald", (1e-10, 1e10), 0, numpy.inf), + ("triangular", (0.0, 1.0, 1.0), 0, 1), + ("vonmises", (-3.0, 1e-3), -numpy.pi, numpy.pi), + ], + ) + def test_extreme_range(self, name, params, lo, hi, dtype): + res = getattr(default_rng(1), name)(*params, size=10**5, dtype=dtype) + res = dpnp.asnumpy(res) + assert ((res >= lo) & (res <= hi)).all() + + @pytest.mark.parametrize( + "name, params", + [ + ("wald", (1.0, 1.0)), + ("pareto", (1.0,)), + ("weibull", (1.0,)), + ("power", (1.0,)), + ("rayleigh", (1.0,)), + ("laplace", (0.0, 1.0)), + ("gumbel", (0.0, 1.0)), + ("logistic", (0.0, 1.0)), + ("vonmises", (0.0, 1.0)), + ("chisquare", (1.0,)), + ("standard_t", (1.0,)), + ("f", (1.0, 1.0)), + ("gamma", (1.0, 1.0)), + ("beta", (1.0, 1.0)), + ], + ) + def test_nan_param(self, name, params): + # NaN passes NumPy's checks and gives NaN; every parameter is tried + rng = default_rng(0) + for i in range(len(params)): + args = list(params) + args[i] = numpy.nan + res = dpnp.asnumpy(getattr(rng, name)(*args, size=100)) + assert numpy.isnan(res).all() + + @pytest.mark.parametrize("dtype", float_dtypes) + @pytest.mark.parametrize( + "a, b", [(1e-3, 1e-3), (0.01, 0.01), (0.3, 1e-4), (50.0, 0.01)] + ) + def test_beta_range(self, a, b, dtype): + # the samples underflow, which devices flushing subnormals to zero + # and with inexact division once turned into inf or 1 + ulp + res = dpnp.asnumpy(default_rng(1).beta(a, b, size=10**5, dtype=dtype)) + assert ((res >= 0) & (res <= 1)).all() + def test_array_params(self): rng = default_rng(3) loc = dpnp.arange(4.0)[:, None] * 10 @@ -255,11 +398,20 @@ def test_array_params(self): dpnp.asnumpy(res.mean(axis=1)), [0.0, 10.0, 20.0, 30.0], atol=0.1 ) + # the output shape is the broadcast shape of the parameters + res = rng.gamma(dpnp.ones(3), numpy.ones((2, 1))) + assert res.shape == (2, 3) + # parameters given as lists res = rng.uniform([0.0, 10.0], [1.0, 11.0], size=(1000, 2)) res = dpnp.asnumpy(res) assert (res[:, 0] < 1).all() and (res[:, 1] >= 10).all() + def test_array_params_chisquare(self): + df = dpnp.asarray([1.0, 10.0]) + res = default_rng(4).chisquare(df, size=(10**5, 2)) + assert_allclose(dpnp.asnumpy(res.mean(axis=0)), [1.0, 10.0], rtol=0.05) + @pytest.mark.parametrize( "size", [ @@ -303,6 +455,13 @@ def test_size_not_broadcastable(self): "name, params, msg", [ ("normal", (0.0, -1.0), "scale < 0"), + ("lognormal", (0.0, -1.0), "sigma < 0"), + ("exponential", (-1.0,), "scale < 0"), + ("standard_gamma", (-1.0,), "shape < 0"), + ("gamma", (1.0, -1.0), "scale < 0"), + ("chisquare", (0.0,), "df <= 0"), + ("beta", (0.0, 1.0), "a <= 0"), + ("beta", (1.0, 0.0), "b <= 0"), ], ) def test_scalar_checks(self, name, params, msg): @@ -314,12 +473,19 @@ def test_uniform_overflow(self): big = numpy.finfo(numpy.float64).max assert_raises(OverflowError, rng.uniform, -big, big) + def test_standard_exponential_method(self): + rng = default_rng(0) + rng.standard_exponential(3, method="inv") + assert_raises(ValueError, rng.standard_exponential, 3, method="x") + @pytest.mark.parametrize("dtype", float_dtypes) @pytest.mark.parametrize( "name, params", [ ("random", ()), ("standard_normal", ()), + ("standard_exponential", ()), + ("standard_gamma", (2.0,)), ], ) def test_out(self, name, params, dtype): @@ -351,6 +517,33 @@ def test_out_errors(self): assert_raises(ValueError, rng.random, out=big[:, ::2]) assert_raises(ValueError, rng.random, out=big[:, :3].T) + def test_device_params_terminate(self): + # the parameters on which the rejection samplers could loop forever + # give NaN + rng = default_rng(0) + nan, inf = numpy.nan, numpy.inf + + res = rng.beta(dpnp.asarray([-inf, -1.0, 0.0, nan]), 1.0, size=(8, 4)) + assert dpnp.isnan(res).all() + + # a shape just below zero would reject nearly every draw + shape = dpnp.asarray([-1e-8, -1e-12, -1e-30, -1.0]) + res = rng.standard_gamma(shape, size=(8, 4)) + assert dpnp.isnan(res).all() + res = rng.chisquare(shape, size=(8, 4)) + assert dpnp.isnan(res).all() + res = rng.standard_gamma(shape.astype(dpnp.float32), dtype="f4") + assert dpnp.isnan(res).all() + rng.standard_t(shape) + + def test_host_params_terminate(self): + # valid parameters on which rounding in single precision made every + # draw of the rejection sampler fail + rng = default_rng(0) + kappa = [5e-5, 1e-4, 3e-4, 4e-4] + res = rng.vonmises(0.0, kappa, size=(8, 4), dtype=dpnp.float32) + assert (dpnp.abs(res) <= numpy.pi).all() + def _assert_same_error(np_call, dp_call): """Assert `dp_call` raises the same exception as `np_call`.""" @@ -377,8 +570,40 @@ class TestChecks: ("normal", (0.0, -1.0)), ("normal", (0.0, -0.0)), ("normal", (0.0, [1.0, -1.0])), + ("lognormal", (0.0, -1.0)), + ("exponential", (-1.0,)), + ("exponential", ([-1.0],)), + ("standard_gamma", (-0.0,)), + ("gamma", (1.0, -1.0)), + ("gamma", ([-1.0], 1.0)), + ("chisquare", (0.0,)), + ("chisquare", ([1.0, 0.0],)), + ("standard_t", (0.0,)), + ("f", (0.0, 1.0)), + ("f", (1.0, [0.0])), + ("beta", (0.0, 1.0)), + ("beta", (1.0, [0.0])), + ("pareto", (0.0,)), + ("power", (0.0,)), + ("weibull", (-1.0,)), + ("weibull", (-0.0,)), + ("rayleigh", (-1.0,)), + ("laplace", (0.0, -1.0)), + ("gumbel", (0.0, -1.0)), + ("logistic", (0.0, -1.0)), + ("wald", (0.0, 1.0)), + ("wald", (1.0, 0.0)), + ("wald", ([1.0], [0.0])), + ("vonmises", (0.0, -1.0)), + ("triangular", (1.0, 0.0, 2.0)), + ("triangular", (0.0, 2.0, 1.0)), + ("triangular", (1.0, 1.0, 1.0)), + ("triangular", ([0.0, 1.0], 1.0, [2.0, 1.0])), ("normal", (0.0, None)), ("uniform", (None, 1.0)), + ("triangular", (0.0, None, 2.0)), + ("exponential", ("x",)), + ("exponential", (1j,)), ], ) @pytest.mark.filterwarnings("ignore:overflow encountered:RuntimeWarning") @@ -394,6 +619,9 @@ def test_same_as_numpy(self, name, args): "name, args", [ ("normal", (None, [1.0, 2.0])), + ("laplace", (0.0, [None, 1.0])), + ("gamma", (None, [1.0, 2.0])), + ("triangular", (None, [1.0, 2.0], 3.0)), ], ) def test_none_with_array_is_nan(self, name, args): @@ -410,6 +638,7 @@ def _draw(self, device, bitgen): rng.random(1000, dtype=dpnp.float32), rng.bit_generator.random_raw(1000), rng.normal(size=1000, dtype=dpnp.float32), + rng.gamma(0.7, size=1000, dtype=dpnp.float32), ) @pytest.mark.parametrize("bitgen", [Philox4x3210]) diff --git a/dpnp/tests/test_sycl_queue.py b/dpnp/tests/test_sycl_queue.py index 4145db0b3dc..b7601cf06a5 100644 --- a/dpnp/tests/test_sycl_queue.py +++ b/dpnp/tests/test_sycl_queue.py @@ -926,10 +926,30 @@ def test_random_state(func, args, kwargs, device, usm_type): _generator_methods = [ + pytest.param("beta", [0.5, 2.0], {"size": (3, 4)}), + pytest.param("chisquare", [3.0], {"size": (3, 4)}), + pytest.param("exponential", [2.0], {"size": (3, 4)}), + pytest.param("f", [2.0, 3.0], {"size": (3, 4)}), + pytest.param("gamma", [2.0], {"size": (3, 4)}), + pytest.param("gumbel", [], {"size": (3, 4)}), + pytest.param("laplace", [], {"size": (3, 4)}), + pytest.param("logistic", [], {"size": (3, 4)}), + pytest.param("lognormal", [], {"size": (3, 4)}), pytest.param("normal", [], {"size": (3, 4)}), + pytest.param("pareto", [2.0], {"size": (3, 4)}), + pytest.param("power", [2.0], {"size": (3, 4)}), pytest.param("random", [], {"size": (3, 4)}), + pytest.param("rayleigh", [], {"size": (3, 4)}), + pytest.param("standard_cauchy", [], {"size": (3, 4)}), + pytest.param("standard_exponential", [], {"size": (3, 4)}), + pytest.param("standard_gamma", [2.0], {"size": (3, 4)}), pytest.param("standard_normal", [], {"size": (3, 4)}), + pytest.param("standard_t", [3.0], {"size": (3, 4)}), + pytest.param("triangular", [0.0, 1.0, 3.0], {"size": (3, 4)}), pytest.param("uniform", [], {"size": (3, 4)}), + pytest.param("vonmises", [0.0, 1.0], {"size": (3, 4)}), + pytest.param("wald", [1.0, 2.0], {"size": (3, 4)}), + pytest.param("weibull", [2.0], {"size": (3, 4)}), ] @@ -964,6 +984,8 @@ def test_bit_generator(bit_generator, device): # queue than the generator's _diff_queue_methods = [ pytest.param("normal", ["f"]), + pytest.param("standard_gamma", ["f"]), + pytest.param("triangular", [0.0, "f", 5.0]), pytest.param("uniform", ["f", 5.0]), ] diff --git a/dpnp/tests/test_usm_type.py b/dpnp/tests/test_usm_type.py index 06e273194e5..8fda5a32b58 100644 --- a/dpnp/tests/test_usm_type.py +++ b/dpnp/tests/test_usm_type.py @@ -1839,6 +1839,8 @@ def test_tensorsolve(self, usm_type, usm_type_other): [ pytest.param("normal", [[1.0, 2.0]], {}), pytest.param("random", [3], {}), + pytest.param("standard_gamma", [2.0], {"size": 3}), + pytest.param("triangular", [0.0, 1.0, 3.0], {"size": 3}), pytest.param("uniform", [], {"size": 3}), ], )