From 66e7a74f9eada7a74cfef3017e6a3158d316c590 Mon Sep 17 00:00:00 2001 From: Nikita Grigorian Date: Wed, 4 Nov 2026 11:46:55 -0800 Subject: [PATCH] Add integers and the discrete distributions to dpnp.random.Generator Add integers and the poisson, binomial, negative_binomial, geometric, logseries, zipf and hypergeometric samplers, along with noncentral_chisquare and noncentral_f, which need poisson. The discrete samplers return int64 and compute in double precision, or in single precision on devices without fp64, where host parameters are clipped to stay inside their domain once rounded. --- CHANGELOG.md | 1 + doc/reference/random.rst | 31 +- dpnp/AGENTS.md | 2 +- dpnp/libtensor/AGENTS.md | 2 +- dpnp/libtensor/extensions/rng/CMakeLists.txt | 5 + dpnp/libtensor/extensions/rng/binomial.cpp | 44 + dpnp/libtensor/extensions/rng/common.hpp | 64 +- dpnp/libtensor/extensions/rng/discrete.cpp | 62 ++ .../extensions/rng/distributions.hpp | 5 + dpnp/libtensor/extensions/rng/integers.cpp | 45 + .../extensions/rng/kernels/distributions.hpp | 795 ++++++++++++++++++ .../extensions/rng/kernels/samplers.hpp | 220 +++++ .../extensions/rng/kernels/umul_hi.hpp | 55 ++ dpnp/libtensor/extensions/rng/noncentral.cpp | 52 ++ dpnp/libtensor/extensions/rng/poisson.cpp | 50 ++ dpnp/libtensor/extensions/rng/rng_py.cpp | 5 + dpnp/libtensor/extensions/rng/sample.cpp | 8 + dpnp/random/_bit_generator.py | 4 +- dpnp/random/_generator.py | 572 ++++++++++++- dpnp/tests/AGENTS.md | 2 +- dpnp/tests/test_random_generator.py | 539 +++++++++++- dpnp/tests/test_sycl_queue.py | 18 + dpnp/tests/test_usm_type.py | 3 + .../cupy/random_tests/test_generator_api.py | 9 - 24 files changed, 2549 insertions(+), 44 deletions(-) create mode 100644 dpnp/libtensor/extensions/rng/binomial.cpp create mode 100644 dpnp/libtensor/extensions/rng/discrete.cpp create mode 100644 dpnp/libtensor/extensions/rng/integers.cpp create mode 100644 dpnp/libtensor/extensions/rng/kernels/umul_hi.hpp create mode 100644 dpnp/libtensor/extensions/rng/noncentral.cpp create mode 100644 dpnp/libtensor/extensions/rng/poisson.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index c2820a39265..15ffd8e3faa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * 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) +* Added `integers` and the discrete distributions of `numpy.random.Generator`, such as `poisson`, `binomial` and `hypergeometric`, to `dpnp.random.Generator`, along with `noncentral_chisquare` and `noncentral_f` [#3159](https://github.com/IntelPython/dpnp/pull/3159) ### Changed diff --git a/doc/reference/random.rst b/doc/reference/random.rst index a44b00f9ce3..f966cb195b9 100644 --- a/doc/reference/random.rst +++ b/doc/reference/random.rst @@ -17,7 +17,7 @@ host-side oneMKL RNG. rng = dpnp.random.default_rng(seed=42) x = rng.standard_normal((3, 4)) # on the default device - y = rng.random(5, usm_type="shared") + y = rng.integers(0, 10, size=5, usm_type="shared") rng = dpnp.random.default_rng(42, device="cpu") @@ -32,7 +32,9 @@ Differences from NumPy type, the raw bits and the uniform floats are the same on every device. The other distributions use transcendental functions, which may round differently on another device: their samples agree to a few ULP, or - differ when a rejection sampler then takes another candidate. + differ when a rejection sampler then takes another candidate. The + discrete distributions also depend on the device's support of double + precision, as noted below. * The values depend on the output size: a draw is not a prefix of a larger one, and two draws of ``n`` values differ from one draw of ``2 * n``. Code splitting its draws into chunks gets other values. @@ -48,7 +50,20 @@ Differences from NumPy not, since that would wait for the device. 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. + memory out of bounds. Floating-point arrays given for integer parameters + are converted as by :meth:`dpnp.ndarray.astype`, whose result for NaN and + out-of-range values is undefined. +* The discrete distributions return ``int64`` values. They compute in double + precision on devices supporting it, and in single precision otherwise. + There, the streams differ from the ones of devices with double precision, + and some samples are only multiples of a power of two: the ones above + ``2**24`` of ``geometric``, ``logseries`` and ``zipf``, and the ones of + ``poisson`` and ``binomial`` with a standard deviation above about + ``2**22``. +* ``zipf`` follows the Zipf distribution truncated to the largest ``int64``. + NumPy's rejection test rounds to accepting every sample above ``2**53``, + which makes its tail too heavy for ``a`` close to 1. +* ``hypergeometric`` requires ``ngood`` and ``nbad`` below ``10**9``. * Generators and bit generators cannot be pickled or copied, and the bit generators have no ``advance`` or ``jumped`` method. @@ -79,6 +94,7 @@ Simple random data :toctree: generated/ :nosignatures: + Generator.integers Generator.random Distributions @@ -89,16 +105,24 @@ Distributions :nosignatures: Generator.beta + Generator.binomial Generator.chisquare Generator.exponential Generator.f Generator.gamma + Generator.geometric Generator.gumbel + Generator.hypergeometric Generator.laplace Generator.logistic Generator.lognormal + Generator.logseries + Generator.negative_binomial + Generator.noncentral_chisquare + Generator.noncentral_f Generator.normal Generator.pareto + Generator.poisson Generator.power Generator.rayleigh Generator.standard_cauchy @@ -111,6 +135,7 @@ Distributions Generator.vonmises Generator.wald Generator.weibull + Generator.zipf Bit generators -------------- diff --git a/dpnp/AGENTS.md b/dpnp/AGENTS.md index 99c14e0b8c0..6009def4f08 100644 --- a/dpnp/AGENTS.md +++ b/dpnp/AGENTS.md @@ -21,7 +21,7 @@ - Extension calls are ordered with `_manager = dpctl.utils.SequentialOrderManager[exec_q]`: pass `depends=_manager.submitted_events`, then `_manager.add_event_pair(ht_ev, ev)`. Don't add blocking waits. - Use `dpnp.default_float_type(...)` or `map_dtype_to_device(...)` for the default float type. - Deprecated ufuncs use `DPNPDeprecatedUnaryFunc(..., deprecated_msg=...)`. -- A `dpnp.random` distribution is a `Generator` method: it calls `self._check_queues(...)` on its array parameters (raises `ExecutionPlacementError`), checks host parameters as NumPy does, and draws through `_cont` or `_float_sample` into `BitGenerator._fill`. Parameters on the device are not checked, as the check waits for the device. With `usm_type=None`, the output and the temporaries take the coerced USM type of the array parameters (`Generator._sample`). +- A `dpnp.random` distribution is a `Generator` method: it calls `self._check_queues(...)` on its array parameters (raises `ExecutionPlacementError`), checks host parameters as NumPy does, and draws through `_cont`, `_float_sample` or `_discrete_sample` into `BitGenerator._fill`. Parameters on the device are not checked, as the check waits for the device. With `usm_type=None`, the output and the temporaries take the coerced USM type of the array parameters (`Generator._sample`). - Unsupported NumPy keywords raise `NotImplementedError` (`dpnp.check_limitations`) and are listed under `Limitations` in the docstring. - The Array API spec version is set in three places that change together: `__array_api_version__` in `dpnp/tensor/_array_api.py`, `ARRAY_API_TESTS_VERSION` in `.github/workflows/conda-package.yml` and `doc/reference/array_api.rst` (plus the version test in `dpnp/tests/test_ndarray.py`). Inspection utilities come from `__array_namespace_info__` in `dpnp_array_api_info.py`. - `dpnp.tensor` replaces `dpctl.tensor`, which no longer exists here; `SyclQueue`, `SyclDevice`, `dpctl.memory` and `SequentialOrderManager` still come from `dpctl`. diff --git a/dpnp/libtensor/AGENTS.md b/dpnp/libtensor/AGENTS.md index f0ca1eb2213..b2e338bba6f 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` draws 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`): raw bits and uniform floats are bit-identical across devices, samplers using transcendental functions agree to a few ULP, or differ where a rejection sampler then takes another candidate. 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. +- `rng` draws 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`): raw bits and uniform floats are bit-identical across devices, samplers using transcendental functions agree to a few ULP, or differ where a rejection sampler then takes another candidate. 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 or the smallest `int64` for the inputs that would make it loop forever, including values that are valid but round badly in fp32. Use `umul_hi` from `kernels/umul_hi.hpp` instead of `sycl::mul_hi`, which DPC++ 2026.2 miscompiles on Intel GPUs. - 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 5faa8f59941..ad41b4a34ac 100644 --- a/dpnp/libtensor/extensions/rng/CMakeLists.txt +++ b/dpnp/libtensor/extensions/rng/CMakeLists.txt @@ -29,9 +29,14 @@ set(python_module_name _rng_impl) set(_module_src + ${CMAKE_CURRENT_SOURCE_DIR}/binomial.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/discrete.cpp ${CMAKE_CURRENT_SOURCE_DIR}/exponential.cpp ${CMAKE_CURRENT_SOURCE_DIR}/gamma.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/integers.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/noncentral.cpp ${CMAKE_CURRENT_SOURCE_DIR}/normal.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/poisson.cpp ${CMAKE_CURRENT_SOURCE_DIR}/rng_py.cpp ${CMAKE_CURRENT_SOURCE_DIR}/sample.cpp ${CMAKE_CURRENT_SOURCE_DIR}/uniform.cpp diff --git a/dpnp/libtensor/extensions/rng/binomial.cpp b/dpnp/libtensor/extensions/rng/binomial.cpp new file mode 100644 index 00000000000..58c8bce687e --- /dev/null +++ b/dpnp/libtensor/extensions/rng/binomial.cpp @@ -0,0 +1,44 @@ +//***************************************************************************** +// 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 discrete_dispatch_table_t binomial_dispatch_table; + +void init_binomial(py::module_ m) +{ + def_discrete_sampler( + m, "_binomial", + "Fill the result with binomial samples of n trials of probability p", + binomial_dispatch_table); +} +} // namespace dpnp::extensions::rng diff --git a/dpnp/libtensor/extensions/rng/common.hpp b/dpnp/libtensor/extensions/rng/common.hpp index 1e72cd95017..abdd80042eb 100644 --- a/dpnp/libtensor/extensions/rng/common.hpp +++ b/dpnp/libtensor/extensions/rng/common.hpp @@ -81,6 +81,10 @@ typedef sycl::event (*sample_fn_ptr_t)(sycl::queue &, using sample_dispatch_table_t = sample_fn_ptr_t[rng_ns::num_engines][td_ns::num_types]; +// Discrete samplers, whose result is always int64, indexed by the engine and +// by whether they compute in double (1) or in float (0) +using discrete_dispatch_table_t = sample_fn_ptr_t[rng_ns::num_engines][2]; + namespace detail { template @@ -138,8 +142,11 @@ sycl::event submit_sample(sycl::queue &exec_q, } // namespace detail // Checks the parameters and the result against SamplerT, then submits the -// kernel, unless there is nothing to sample. -template +// kernel, unless there is nothing to sample. ComputeT is the type the +// sampler computes in, which the device must support. +template sycl::event sample_impl(sycl::queue &exec_q, const rng_ns::EngineState &state, const std::size_t nelems, @@ -162,6 +169,7 @@ sycl::event sample_impl(sycl::queue &exec_q, } detail::validate_type_for_device(exec_q); + detail::validate_type_for_device(exec_q); detail::validate_types_for_device( exec_q, std::make_index_sequence{}); @@ -195,6 +203,12 @@ struct RealTypes { }; +// Types the integer samplers support +template +struct IntegerTypes : std::bool_constant> +{ +}; + // The raw output of a bit generator, whole or as 32-bit halves template struct RawTypes : std::bool_constant || @@ -238,11 +252,27 @@ void init_sample_dispatch_table(sample_dispatch_table_t &table) table[philox]); } +template