Repository navigation
Add dpnp.random.Generator on the oneMKL device RNG API - #3157
ndgrigorian wants to merge 11 commits into
Conversation
dbdb2c8 to
3963d20
Compare
3963d20 to
6c3ff46
Compare
|
View rendered docs @ https://intelpython.github.io/dpnp/pull/3157/index.html |
|
Array API standard conformance tests for dpnp=0.22.0dev1=np2py314h8d9cdd5_16 ran successfully. |
There was a problem hiding this comment.
🟡 Changes recommended
Unresolved RNG correctness issues can deadlock state access, duplicate random streams through races or invalid sizes, and produce incorrect results for large draws or overflowing float32 ranges.
5 open issues
What changed in this PR
Introduces the foundation for device-side random generation through NumPy’s Generator API, the first step in the RNG migration series. Legacy RandomState and module-level functions remain unchanged.
Changes:
- Adds
Generator,default_rng, and Philox bit generators with state management and spawning. - Implements base distributions in a header-only oneMKL device RNG extension, with strided and GPU-optimized blocked sampling.
- Adds API documentation and tests for distributions, reproducibility, queues, and allocation types.
| File | Description |
|---|---|
dpnp/tests/test_usm_type.py |
Tests RNG allocation types. |
dpnp/tests/test_sycl_queue.py |
Tests queue placement and mismatch errors. |
dpnp/tests/test_random_generator.py |
Tests generators, distributions, state, and reproducibility. |
dpnp/tests/AGENTS.md |
Adds RNG testing guidance. |
dpnp/random/_generator.py |
Implements Generator and base distributions. |
dpnp/random/_bit_generator.py |
Implements Philox state, raw draws, and spawning. |
dpnp/random/__init__.py |
Exports the new random APIs. |
dpnp/libtensor/extensions/rng/uniform.cpp |
Registers uniform and raw-bit samplers. |
dpnp/libtensor/extensions/rng/sample.cpp |
Validates inputs and manages sampling lifetimes. |
dpnp/libtensor/extensions/rng/rng_py.cpp |
Defines the Python extension module. |
dpnp/libtensor/extensions/rng/normal.cpp |
Registers the normal sampler. |
dpnp/libtensor/extensions/rng/kernels/samplers.hpp |
Defines distribution sampler adapters. |
dpnp/libtensor/extensions/rng/kernels/sample.hpp |
Implements strided and blocked sampling layouts. |
dpnp/libtensor/extensions/rng/kernels/engines.hpp |
Wraps Philox engines and uniform-bit conversion. |
dpnp/libtensor/extensions/rng/kernels/distributions.hpp |
Implements parameters and Box–Muller sampling. |
dpnp/libtensor/extensions/rng/distributions.hpp |
Declares sampler registration functions. |
dpnp/libtensor/extensions/rng/common.hpp |
Provides dispatch and kernel submission helpers. |
dpnp/libtensor/extensions/rng/CMakeLists.txt |
Configures the header-only RNG extension. |
dpnp/libtensor/AGENTS.md |
Documents RNG architecture and stream invariants. |
dpnp/CMakeLists.txt |
Includes the RNG extension in the build. |
dpnp/AGENTS.md |
Adds Python RNG implementation guidance. |
doc/reference/random.rst |
Documents modern RNG APIs alongside legacy functions. |
CHANGELOG.md |
Records the new device RNG API. |
Review effort: Balanced
This comment was marked as off-topic.
This comment was marked as off-topic.
antonwolfy
left a comment
There was a problem hiding this comment.
The sampling core is right: the stream matches Philox4x32-10 bit for bit on GPU and CPU, and it passes PractRand and TestU01. The problems are at the edges, in the Python layer, the docs and the tests. Copilot's five points all hold; of the inline comments below, I'd fix the overlap check, the device claim in the docs and the skipped CuPy tests before merge.
Open bot points I confirm:
- Copilot,
lockdeadlock (_bit_generator.py:110). New evidence: a thread runningwith bg.lock: bg.statewas still alive after a 3 s join. - Copilot,
random_raw(output=False)rewind (_bit_generator.py:242). New evidence:random_raw(16), thenrandom_raw(-16, output=False), thenrandom_raw(16)returns the same 16 values, while NumPy'sPhilox.random_raw(-16, output=False)leaves the counter unchanged.random_raw((2, -8), output=False)rewinds too, andrandom_raw((2**40, 2**40), output=False)wrapsnumpy.prodto 0 and leaves the state alone. - Copilot, key read after the lock is released (
_bit_generator.py:391). - Copilot,
LARGE_RANGES(rng/CMakeLists.txt:40). It takes more than 1.7e10 elements to reach. - Copilot, float32
uniformrange (_generator.py:129). New evidence:uniform(-3e38, 3e38, size=1000, dtype=numpy.float32)returns 1000infs.
What I ran, on a Panther Lake U iGPU (Level Zero and OpenCL) and the OpenCL CPU device, with the PR built in place:
test_random_generator.py: 102 passed on level_zero:gpu; 100 passed and 2 skipped (the two-device tests) on opencl:cpu. The newtest_sycl_queue.pyandtest_usm_type.pytests,third_party/cupy/random_tests/as shipped,test_random.pyandtest_random_state.pypass on both, with no pytest warnings. The build had no compiler warnings.- A NumPy implementation of Philox4x32-10, checked against the Random123 known-answer vectors, reproduces
random_rawbit for bit in 274 checks on both devices: sizes from 1 to 220+3, a counter carrying into the high word, the wrap at 2128, consecutive calls andoutput=False.randomanduniformmatch its float conversion bit for bit. - PractRand 0.94: 256 GB of the GPU stream with no anomaly. 4-16 GB each, also with none: the CPU stream; keys 0-3 with the same counter, seeds 0-3 and four
spawnchildren, each interleaved by word and by chunk; and streams of 1-, 9- and 17-value calls (built with the reference, which matches the device). - TestU01 1.2.3: SmallCrush passes on the GPU stream, the CPU stream and the interleaved related keys. Crush passes all 144 tests on two seeds; a third seed had p = 6.4e-4 (CouponCollector) and 0.9992 (CollisionPermut), which the other two seeds didn't repeat. BigCrush is still running.
- Distributions, 1e8 draws per case on both devices: KS and chi-square for
random,uniform(ranges up to ±1e30 in float32 and ±1e300 in float64),standard_normalandnormal, in float32 and float64, all in range. Over 6.7e9 normals there is no correlation between the two values of a Box-Muller pair or between neighbouring work-items.
Not inline:
uniformskips the checks of its host parameters when any parameter is on the device.uniform(dpnp.zeros(1), numpy.array([numpy.inf]))returns[inf], where NumPy raisesOverflowError: Range exceeds valid bounds._check_arrays(_generator.py:142-152) returnsNonefor every parameter as soon as one is a device array, whilenormalchecks each parameter on its own (normal(dpnp.zeros(1), -1.0)raises). Suggestion: check that the host parameters are finite even when the range can't be computed on the host.- The output ignores the parameters' USM type, and the new
dpnp/AGENTS.mdline says something else about temporaries.normal(dpnp.zeros(4, usm_type="shared"), 1.0).usm_typeis"device". The line says "Temporaries take the output'susm_type", but_to_kernel_param(_generator.py:360) callsdpnp.asarraywithoutusm_type, so a parameter given as a list becomes a"device"array even withusm_type="shared". Suggestion: passusm_type=usm_typethere, and either follow the inputs' USM type as other dpnp functions do or list the keyword-only rule under "Differences from NumPy". - The method docstrings don't describe
dtype,usm_typeandout. They appear only as one-lineLimitations, and the methods have noParameters,ReturnsorExamples, unlike the legacydpnp.random.uniform.random.rstdoesn't listGenerator.spawn,Generator.bit_generator,BitGenerator.random_raworBitGenerator.state, and thedevicetext at_bit_generator.py:87refers todpctl.tensor.Device, which isdpnp.tensor.Deviceelsewhere in dpnp now. - Nits:
- Nit: the
_random_rawdoc string (uniform.cpp:46-47) says "low half first", butBitGen::next_uint64returns(first_word << 32) | second_word, so the first 32-bit word is the high half. - Nit:
NormalSampler::operator()(kernels/samplers.hpp) lacks the#pragma clang fp contract(off)thatUniformSamplerhas, soloc + scale * zmay become an FMA on one device and not on another (icpx defaults to-ffp-contract=fast-honor-pragmas). - Nit: the Sphinx spelling step reports 41 new words from these docstrings (Philox, BitGenerator, SeedSequence, ziggurat, ...); they belong in
doc/known_words.txt. - Non-blocking:
pickleandcopy.deepcopyof aGeneratoror aBitGeneratorfail withTypeError(from theSyclQueueand theLock), andadvance/jumpedare absent; NumPy has all of them. Worth a line under "Differences from NumPy" until a follow-up adds them.
- Nit: the
Checked, no issues
- Engine setup:
{key, counter_lo, counter_hi + w}maps to oneMKL's key k0/k1 and counter c0-c3 the way_launchassumes, and onegenerateonphilox4x32x10<8>or<16>yields the words of 8 or 16 calls on<1>. Raw and uniform values are bit-identical between the GPU (blocked) and CPU (strided) layouts. - Stream separation: each work-item owns 264 counters, and the carry from the low word keeps the blocks disjoint up to the wrap at 2128.
- Lifetimes and validation:
keep_args_aliveholds the result and the parameters until the kernel ends; the queue, writability, contiguity, size, dtype, engine and seed-length checks are in place. - fp64: the double kernels are guarded and split per kernel, and the float paths use only float literals.
- NumPy parity: size normalisation, dtype handling, the
scale < 0,-0.0, NaN andNonechecks, broadcasting and theExecutionPlacementErrorcases match NumPy. - Statistical tests: the mean and variance checks of
test_momentsgave no failure over 200 other seeds. - CI: GitHub Actions and the required Jenkins lanes are green. The failed Jenkins Windows py3.14 wheel build is an icx out-of-memory crash in
copy_and_cast_usm_to_usm.cpp, which this PR doesn't touch.
9fef498 to
1bd3cee
Compare
|
Thanks for the thorough review. All of the inline points are addressed in cab62a2, with replies on each thread. On the points from the review body:
The six PRs above are rebased onto this one, carrying the changes where they apply: |
|
BigCrush (TestU01 1.2.3) has finished on the GPU stream: all 160 statistics of its 106 tests passed, over 3.57e11 32-bit words. That completes the randomness testing from my review: no battery found a problem with the generator. |
antonwolfy
left a comment
There was a problem hiding this comment.
Re-checked at this head on level_zero:gpu and opencl:cpu: all the earlier points are fixed, and the tests pass on both devices. A few new points inline: the first two are worth fixing, the rest are nits.
There was a problem hiding this comment.
🟡 Changes recommended
Host-parameter uploads still block repeated sampling, and the enabled uniform KS test fails on devices without fp64 because the NumPy reference dtype does not match the sampling dtype.
2 open issues
5 resolved since last review
Review effort: Balanced
antonwolfy
left a comment
There was a problem hiding this comment.
Re-checked at this head on a Panther Lake H (Arc B390, fp64), on level_zero:gpu and opencl:cpu, against master at the merge-base. All the points of the last round are fixed, and the oneMath switch looks right. Both of Copilot's new points hold; I replied on those threads with measurements. A few small points inline.
Tests. The full suite, one pytest process per file, passes on both devices: 94505 passed on level_zero:gpu and 94525 on opencl:cpu, with no failures or crashes. It gives no new warnings: the only one, invalid value encountered in sign from CuPy's test_sign_nan, comes on master too. With both devices visible, test_random_generator.py passes 138 of 138, and the RNG test files pass 20 repeats on both devices.
Build and footprint (clean build, 16 cores, default ninja jobs):
| master | this PR | |
|---|---|---|
| Wall time | 835 s | 837 s |
Peak RSS of one process (ld.lld) |
3.40 GiB | 3.40 GiB |
| Peak RSS of all processes | 9.26 GiB | 9.31 GiB |
| Compiler and CMake warnings, legacy backend included | 0 | 0 |
| Extension modules on disk | 183.4 MiB | 184.3 MiB |
import dpnp |
347 ms, 405 MB RSS | 355 ms, 416 MB RSS |
_rng_impl takes 59 CPU-seconds, is done long before the _tensor_impl link that ends the build, and links no oneMKL library. Most of the import cost is numpy.random and the libcrypto it loads, which _bit_generator.py now imports at import time for SeedSequence.
Performance against the legacy RandomState (MCG59 on the GPU, MT19937 on the CPU), median ms per call for 1e7 samples:
| device | method | Generator float32 | Generator float64 | RandomState float64 |
|---|---|---|---|---|
| level_zero:gpu | random |
0.67 | 1.29 | 0.98 |
| level_zero:gpu | standard_normal |
0.86 | 3.14 | 2.95 |
| level_zero:gpu | normal |
0.96 | 3.21 | 2.95 |
| level_zero:gpu | uniform |
0.86 | 1.33 | 0.96 |
| opencl:cpu | random |
6.5 | 8.6 | 26.6 |
| opencl:cpu | standard_normal |
8.6 | 10.5 | 44.6 |
| opencl:cpu | normal |
8.3 | 10.5 | 45.0 |
| opencl:cpu | uniform |
5.3 | 8.5 | 26.1 |
- On the GPU, float64 runs at 0.6-1.0x of
RandomStatefrom 1e5 to 1e7 samples, since Philox costs more than MCG59. On the CPU the Generator is 2.2-7.6x faster from 1e6 samples up. NumPy'sGenerator(Philox)takes 53-132 ms for the same float64 draws. - A small call (1 to 100 samples) takes 35-63 µs on the GPU and 43-95 µs on the CPU, against 15-23 µs for
RandomState. It is all host time, spread over theSequentialOrderManagerbookkeeping, the extension call with its keep-alive host task, and the NumPy shape checks, which run for scalar parameters too. -fno-fast-mathcosts nothing measurable: without it, the GPU times stay within noise and the CPU is 3-9% slower at 1e8.- The legacy API's timings are unchanged.
Coverage: +0.11% (84.397% to 84.508%). _bit_generator.py is at 100% and _generator.py at 98.9%, missing only the C-contiguous out check left for #3158.
Add the dpnp.libtensor._rng_impl extension, which samples on the device from oneMKL's device philox4x32x10 engine, and the Generator API of NumPy on top of it: BitGenerator, Philox4x3210, SeedSequence (NumPy's), default_rng, and the random, uniform, standard_normal and normal methods. A bit generator keeps its state on the host and every launch rebuilds its engines from it, so the values depend only on the seed and the output size, not on the device. The legacy RandomState and module functions are unchanged.
- Make the bit generator lock reentrant and read the Philox key under it - Validate `size` in `random_raw(output=False)` and count values without overflow - Drop the seed sequence of a Philox4x3210 given its key, as NumPy does - Check the `uniform` range in the sampling dtype, and host bounds next to device ones - Raise TypeError for parameters that don't cast to the sampling dtype - Follow the USM type of the array parameters when `usm_type` is None - Accept a usm_ndarray and, for `random` and `standard_normal`, an F-contiguous `out` - Reject a parameter array overlapping the result, and build the module with LARGE_RANGES - Take the float32 Box-Muller u1 from a whole 32-bit word, and keep FMA contraction off in the normal sampler - Correct the device-independence claims, document the output-size dependence, and expand the method docstrings - Enable the CuPy Generator and Philox tests, skipping the classes of later distributions and bit generators
- Build `_rng_impl` with -fno-fast-math, so `uniform` over a tiny range is not reassociated on CPU devices and stays below `high` - Pin the float32 Box-Muller tail with a counter whose u1 word is 6 - Test the `bit_generator` setter error, a dpnp-array seed, `sycl_device`, repr and str, and the base BitGenerator methods - Drop the NumPy notice from the distributions header - Leave the kernel headers of the extensions out of the C++ coverage - Use the `dpnp.empty` text for `device` and `sycl_queue`, describe the `Generator.spawn` parameters, and list the generator properties in the docs - Correct the `random_raw(output=False)`, copying and 32-bit raw output descriptions, and rename the legacy "Random generator" section
The rng extension takes the header-only device RNG API from oneMath, with its include directories, when dpnp is built with oneMath, as the blas, fft, lapack and sparse extensions do.
9f9f13e to
88c286e
Compare


This PR is the first of a series splitting #3147 into smaller pieces. It adds the new framework for device-side random number generation and NumPy's
GeneratorAPI, with only the base distributionsThe new
dpnp.libtensor._rng_implextension samples on the device using the oneMKL device RNG API (header-only, so no oneMKL RNG library is linked). A bit generator keeps its state on the host, and each kernel launch rebuilds its engines from that state, so sampling never waits on the device. For a given seed, output size and dtype, the raw bits and the uniform floats are the same on every device, and the normal samples agree to a few ULPAdds
dpnp.random.Generatorwithrandom,uniform,standard_normal,normalandspawn, anddefault_rng; withusm_type=None(the default) the samples take the coerced USM type of the array parametersdpnp.random.BitGenerator(withrandom_raw,spawn,stateandlock),SeedSequence(NumPy's) and thePhilox4x3210bit generator on oneMKL's devicephilox4x32x10The legacy
RandomStateand the module-level functions are untouched until later in the seriesSeries:
dpnp.random.Generatoron the oneMKL device RNG API #3157 (this PR): framework,Philox4x3210and the base distributionsdpnp.random.Generator#3158: Continuous distributionsdpnp.random.Generator#3159: Discrete distributionsdpnp.random.Generator#3160: Multivariate distributions and permutationsPCG64DXSMandMRG32k3abit generators todpnp.random#3161:PCG64DXSMandMRG32k3abit generatorsdpnp.random.RandomStateon the device RNG #3162:RandomStateand the module functions on the device RNG