diff --git a/.flake8 b/.flake8 index a8e21cf..10fccbd 100644 --- a/.flake8 +++ b/.flake8 @@ -28,6 +28,7 @@ per-file-ignores = mkl_random/interfaces/__init__.py: F401 mkl_random/tests/*.py: D102 mkl_random/tests/**/*.py: D102 + benchmarks/**/*.py: D102 filename = *.py, *.pyx, *.pxi, *.pxd max_line_length = 80 diff --git a/.gitignore b/.gitignore index d7e11b0..bd23b38 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,6 @@ mkl_random/mklrand.cpp # Byte-compiled / optimized / DLL files __pycache__/ + +# ASV benchmark artifacts +.asv/ diff --git a/CHANGELOG.md b/CHANGELOG.md index 5f00d69..ecf2d23 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * Added support for `array_like` (broadcastable) `low`/`high` bounds in `randint` [gh-168](https://github.com/IntelPython/mkl_random/pull/168) * Added support for free-threaded (GIL-disabled) CPython builds: the Cython extension is compiled with `freethreading_compatible=True`, so importing `mkl_random` no longer re-enables the GIL [gh-159](https://github.com/IntelPython/mkl_random/pull/159) * Added a thread-safety section to the how-to guide for free-threaded Python [gh-159](https://github.com/IntelPython/mkl_random/pull/159) +* Added an [ASV](https://asv.readthedocs.io/en/stable/) benchmark suite under `benchmarks/` and a `benchmark` optional dependency group [gh-184](https://github.com/IntelPython/mkl_random/pull/184) ### Changed * Sped up `normal`, `uniform`, `exponential`, `laplace`, `gumbel`, `logistic`, `rayleigh` and `lognormal` for array-valued parameters. Seeded results change for these array paths; scalar paths are unchanged [gh-171](https://github.com/IntelPython/mkl_random/pull/171) diff --git a/benchmarks/README.md b/benchmarks/README.md new file mode 100644 index 0000000..fd76fb9 --- /dev/null +++ b/benchmarks/README.md @@ -0,0 +1,107 @@ +# mkl_random ASV Benchmarks + +Performance benchmarks for [mkl_random](https://github.com/IntelPython/mkl_random) using +[Airspeed Velocity (ASV)](https://asv.readthedocs.io/en/stable/). + +### Coverage + +| File | API | Cases | Engines | Sizes | +|------|-----|-------|---------|-------| +| `bench_continuous.py` | `MKLRandomState` | Every continuous distribution; shape parameters chosen to hit each oneMKL method branch (e.g. gamma a>1, 0.6= _MIN_THREADS else "1" + + +_THREADS = os.environ.get("MKL_NUM_THREADS", _thread_count()) +os.environ["MKL_NUM_THREADS"] = _THREADS diff --git a/benchmarks/benchmarks/_utils.py b/benchmarks/benchmarks/_utils.py new file mode 100644 index 0000000..31bf379 --- /dev/null +++ b/benchmarks/benchmarks/_utils.py @@ -0,0 +1,64 @@ +# Copyright (c) 2026, Intel Corporation +# +# 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 Intel Corporation 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 OWNER 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. + +"""Shared helpers and parameter axes for mkl_random benchmarks.""" + +import mkl_random + +_SEED = 42 + +# size axis shared across multiple files: per-call overhead and throughput +_SIZES = [1_000, 1_000_000] + +# Every deterministic basic generator. NONDETERM reads a hardware entropy +# source, so its timings are neither reproducible nor available everywhere. +_ENGINES = [ + "MT19937", + "SFMT19937", + "MT2203", + "MCG31", + "MCG59", + "MRG32K3A", + "R250", + "WH", + "PHILOX4X32X10", + "ARS5", +] + +# Generators that implement viRngUniformBits32/64, which full-range integer +# fills and bytes() call directly. +_ENGINES_BITS = [ + "MT19937", + "SFMT19937", + "MT2203", + "MCG59", + "PHILOX4X32X10", + "ARS5", +] + + +def _make_state(brng="MT19937"): + """Return a seeded MKLRandomState for the given basic generator.""" + return mkl_random.MKLRandomState(_SEED, brng=brng) diff --git a/benchmarks/benchmarks/bench_array_params.py b/benchmarks/benchmarks/bench_array_params.py new file mode 100644 index 0000000..8a795b8 --- /dev/null +++ b/benchmarks/benchmarks/bench_array_params.py @@ -0,0 +1,81 @@ +# Copyright (c) 2026, Intel Corporation +# +# 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 Intel Corporation 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 OWNER 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. + +"""Benchmarks for array-valued distribution parameters. + +One parameter value per output element takes a separate code path from +scalar parameters, so these are tracked apart from bench_continuous.py and +bench_discrete.py. +""" + +import numpy as np + +from ._utils import _SEED, _make_state + +# Smaller than _utils._SIZES: the per-element path is much slower than a +# scalar fill. +_ARRAY_SIZES = [1_000, 100_000] + + +class ArrayParams: + """Distributions called with one parameter value per output element.""" + + params = [_ARRAY_SIZES] + param_names = ["size"] + + def setup(self, size): + self.rs = _make_state() + rng = np.random.default_rng(_SEED) + self.loc = rng.standard_normal(size) + self.scale = rng.uniform(0.5, 2.0, size) + self.high = self.loc + self.scale + self.shape = rng.uniform(0.5, 5.0, size) + self.lam = rng.uniform(1.0, 100.0, size) + self.n = rng.integers(1, 100, size) + self.p = rng.uniform(0.1, 0.9, size) + self.rs.normal(self.loc, self.scale) + self.rs.uniform(self.loc, self.high) + self.rs.exponential(self.scale) + self.rs.standard_gamma(self.shape) + self.rs.poisson(self.lam) + self.rs.binomial(self.n, self.p) + + def time_normal(self, size): + self.rs.normal(self.loc, self.scale) + + def time_uniform(self, size): + self.rs.uniform(self.loc, self.high) + + def time_exponential(self, size): + self.rs.exponential(self.scale) + + def time_standard_gamma(self, size): + self.rs.standard_gamma(self.shape) + + def time_poisson(self, size): + self.rs.poisson(self.lam) + + def time_binomial(self, size): + self.rs.binomial(self.n, self.p) diff --git a/benchmarks/benchmarks/bench_continuous.py b/benchmarks/benchmarks/bench_continuous.py new file mode 100644 index 0000000..e8d52b8 --- /dev/null +++ b/benchmarks/benchmarks/bench_continuous.py @@ -0,0 +1,103 @@ +# Copyright (c) 2026, Intel Corporation +# +# 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 Intel Corporation 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 OWNER 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. + +"""Benchmarks for continuous distributions of mkl_random.MKLRandomState. + +Each distribution gets its own benchmark class (Bench_) so ASV renders a +separate grid tile per distribution. Params are size only. Arguments are +chosen so that every distinct sampling code path is covered. +""" + +from ._utils import _SIZES, _make_state + +# name -> (MKLRandomState method, positional arguments before ``size``) +_CONTINUOUS = { + "uniform": ("uniform", (0.0, 1.0)), + "random_sample": ("random_sample", ()), + "standard_normal": ("standard_normal", ()), + "normal": ("normal", (1.0, 2.0)), + "standard_exponential": ("standard_exponential", ()), + "exponential": ("exponential", (2.0,)), + # MKL's GNORM gamma method may select different internal algorithms + # around shape 1 and 0.6; these shapes exercise those regimes + "standard_gamma": ("standard_gamma", (3.0,)), + "standard_gamma_mid_shape": ("standard_gamma", (0.8,)), + "standard_gamma_small_shape": ("standard_gamma", (0.5,)), + # MKL's CJA beta method may select different internal algorithms (Cheng, + # Johnk or Atkinson) depending on shapes; these shapes exercise those + # regimes + "beta": ("beta", (2.0, 5.0)), + "beta_small_shape": ("beta", (0.5, 0.5)), + "beta_mixed_shape": ("beta", (0.5, 2.0)), + "chisquare": ("chisquare", (3.0,)), + "standard_t": ("standard_t", (5.0,)), + "standard_cauchy": ("standard_cauchy", ()), + "lognormal": ("lognormal", (0.0, 1.0)), + "laplace": ("laplace", (0.0, 1.0)), + "gumbel": ("gumbel", (0.0, 1.0)), + "logistic": ("logistic", (0.0, 1.0)), + "rayleigh": ("rayleigh", (1.0,)), + "wald": ("wald", (1.0, 2.0)), + "weibull": ("weibull", (1.5,)), + "pareto": ("pareto", (3.0,)), + "power": ("power", (2.0,)), + "triangular": ("triangular", (0.0, 0.5, 1.0)), + # von Mises has separate samplers for kappa <= 1 and kappa > 1 + "vonmises": ("vonmises", (0.0, 4.0)), + "vonmises_small_kappa": ("vonmises", (0.0, 0.5)), + "f": ("f", (5.0, 10.0)), + # noncentral chi-square has separate samplers for df > 1 and df < 1 + "noncentral_chisquare": ("noncentral_chisquare", (3.0, 2.0)), + "noncentral_chisquare_small_df": ("noncentral_chisquare", (0.5, 2.0)), + "noncentral_f": ("noncentral_f", (5.0, 10.0, 2.0)), +} + + +def _make_bench(name, method, args): + def setup(self, size): + self._draw = getattr(_make_state(), method) + self._draw(*args, size=size) + + def time_sample(self, size): + self._draw(*args, size=size) + + return type( + f"Bench_{name}", + (), + { + "params": [_SIZES], + "param_names": ["size"], + "setup": setup, + "time_sample": time_sample, + }, + ) + + +globals().update( + { + f"Bench_{name}": _make_bench(name, method, args) + for name, (method, args) in _CONTINUOUS.items() + } +) diff --git a/benchmarks/benchmarks/bench_discrete.py b/benchmarks/benchmarks/bench_discrete.py new file mode 100644 index 0000000..3a92eda --- /dev/null +++ b/benchmarks/benchmarks/bench_discrete.py @@ -0,0 +1,79 @@ +# Copyright (c) 2026, Intel Corporation +# +# 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 Intel Corporation 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 OWNER 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. + +"""Benchmarks for discrete distributions of mkl_random.MKLRandomState. + +Each distribution gets its own benchmark class (Bench_) so ASV renders a +separate grid tile per distribution. Params are size only. ``randint`` is +covered in bench_integers.py. +""" + +from ._utils import _SIZES, _make_state + +# name -> (MKLRandomState method, positional arguments before ``size``) +_DISCRETE = { + # MKL's BTPE binomial method may select a different internal algorithm + # when n * min(p, 1 - p) >= 30; these parameters exercise both regimes + "binomial": ("binomial", (10, 0.5)), + "binomial_large_n": ("binomial", (1000, 0.3)), + "negative_binomial": ("negative_binomial", (5, 0.5)), + # default method (POISNORM); per-method timings are in bench_methods.py + "poisson": ("poisson", (10.0,)), + "geometric": ("geometric", (0.3,)), + # MKL's H2PE hypergeometric method may select a different internal + # algorithm for a large mode; these parameters exercise both regimes + "hypergeometric": ("hypergeometric", (10, 20, 5)), + "hypergeometric_large_mode": ("hypergeometric", (1000, 2000, 500)), + "zipf": ("zipf", (2.0,)), + "logseries": ("logseries", (0.9,)), +} + + +def _make_bench(name, method, args): + def setup(self, size): + self._draw = getattr(_make_state(), method) + self._draw(*args, size=size) + + def time_sample(self, size): + self._draw(*args, size=size) + + return type( + f"Bench_{name}", + (), + { + "params": [_SIZES], + "param_names": ["size"], + "setup": setup, + "time_sample": time_sample, + }, + ) + + +globals().update( + { + f"Bench_{name}": _make_bench(name, method, args) + for name, (method, args) in _DISCRETE.items() + } +) diff --git a/benchmarks/benchmarks/bench_engines.py b/benchmarks/benchmarks/bench_engines.py new file mode 100644 index 0000000..a6f83f4 --- /dev/null +++ b/benchmarks/benchmarks/bench_engines.py @@ -0,0 +1,76 @@ +# Copyright (c) 2026, Intel Corporation +# +# 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 Intel Corporation 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 OWNER 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. + +"""Benchmarks across the basic random number generators (``brng``).""" + +import numpy as np + +from ._utils import _ENGINES, _ENGINES_BITS, _make_state + +_N = 1_000_000 + + +class EngineFill: + """Uniform, Gaussian and narrow-range integer fills for every engine.""" + + params = [_ENGINES] + param_names = ["brng"] + + def setup(self, brng): + self.rs = _make_state(brng) + self.rs.random_sample(_N) + self.rs.standard_normal(_N) + self.rs.randint(0, 1000, _N, dtype=np.int32) + + def time_random_sample(self, brng): + self.rs.random_sample(_N) + + def time_standard_normal(self, brng): + self.rs.standard_normal(_N) + + def time_randint_int32(self, brng): + self.rs.randint(0, 1000, _N, dtype=np.int32) + + +class EngineBits: + """Raw-bit fills, for the engines that implement viRngUniformBits.""" + + params = [_ENGINES_BITS] + param_names = ["brng"] + + def setup(self, brng): + self.rs = _make_state(brng) + self.rs.randint(0, 2**32, _N, dtype=np.uint32) + self.rs.randint(0, 2**64, _N, dtype=np.uint64) + self.rs.bytes(_N) + + def time_randint_uint32_full(self, brng): + self.rs.randint(0, 2**32, _N, dtype=np.uint32) + + def time_randint_uint64_full(self, brng): + self.rs.randint(0, 2**64, _N, dtype=np.uint64) + + def time_bytes(self, brng): + self.rs.bytes(_N) diff --git a/benchmarks/benchmarks/bench_integers.py b/benchmarks/benchmarks/bench_integers.py new file mode 100644 index 0000000..2b3d9f4 --- /dev/null +++ b/benchmarks/benchmarks/bench_integers.py @@ -0,0 +1,95 @@ +# Copyright (c) 2026, Intel Corporation +# +# 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 Intel Corporation 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 OWNER 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. + +"""Benchmarks for randint across dtypes, ranges and bound shapes.""" + +import numpy as np + +from ._utils import _SEED, _SIZES, _make_state + +# case -> (dtype, low, high). The ranges select the integer fill paths: +# small narrow range, drawn with viRngUniform +# wide non-power-of-two range at or above INT_MAX (shifted 32-bit draw, +# or masked 64-bit draw with rejection) +# pow2 power-of-two range above INT_MAX (masked 64-bit draw, no rejection) +# full whole dtype range (raw bits for 32- and 64-bit dtypes) +_CASES = { + "bool": ("bool", 0, 2), + "int8_small": ("int8", 0, 100), + "int8_full": ("int8", -(2**7), 2**7), + "uint8_small": ("uint8", 0, 100), + "uint8_full": ("uint8", 0, 2**8), + "int16_small": ("int16", 0, 100), + "int16_full": ("int16", -(2**15), 2**15), + "uint16_small": ("uint16", 0, 100), + "uint16_full": ("uint16", 0, 2**16), + "int32_small": ("int32", 0, 100), + "int32_wide": ("int32", -(2**30), 2**31), + "int32_full": ("int32", -(2**31), 2**31), + "uint32_small": ("uint32", 0, 100), + "uint32_wide": ("uint32", 0, 3 * 2**30), + "uint32_full": ("uint32", 0, 2**32), + "int64_small": ("int64", 0, 100), + "int64_wide": ("int64", 0, 3 * 2**32), + "int64_pow2": ("int64", 0, 2**40), + "int64_full": ("int64", -(2**63), 2**63), + "uint64_small": ("uint64", 0, 100), + "uint64_wide": ("uint64", 0, 3 * 2**32), + "uint64_pow2": ("uint64", 0, 2**40), + "uint64_full": ("uint64", 0, 2**64), +} + +_BROADCAST_DTYPES = ["uint8", "int32", "int64"] + + +class Randint: + """randint with scalar bounds.""" + + params = [list(_CASES), _SIZES] + param_names = ["case", "size"] + + def setup(self, case, size): + self.rs = _make_state() + self.dtype, self.low, self.high = _CASES[case] + self.rs.randint(self.low, self.high, size, dtype=self.dtype) + + def time_randint(self, case, size): + self.rs.randint(self.low, self.high, size, dtype=self.dtype) + + +class RandintBroadcast: + """randint with an array of upper bounds (one bound per element).""" + + params = [_BROADCAST_DTYPES, _SIZES] + param_names = ["dtype", "size"] + + def setup(self, dtype, size): + self.rs = _make_state() + rng = np.random.default_rng(_SEED) + self.high = rng.integers(2, 100, size).astype(dtype) + self.rs.randint(0, self.high, dtype=dtype) + + def time_randint(self, dtype, size): + self.rs.randint(0, self.high, dtype=dtype) diff --git a/benchmarks/benchmarks/bench_interfaces.py b/benchmarks/benchmarks/bench_interfaces.py new file mode 100644 index 0000000..17b8b7f --- /dev/null +++ b/benchmarks/benchmarks/bench_interfaces.py @@ -0,0 +1,103 @@ +# Copyright (c) 2026, Intel Corporation +# +# 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 Intel Corporation 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 OWNER 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. + +"""Benchmarks for the NumPy drop-in paths. + +``mkl_random.interfaces.numpy_random`` module functions, and ``numpy.random`` +while patched by ``mkl_random.patch_numpy_random()``. +""" + +import numpy as np + +import mkl_random +from mkl_random.interfaces import numpy_random + +from ._utils import _SEED, _SIZES + + +class NumpyRandomInterface: + """Module functions of mkl_random.interfaces.numpy_random.""" + + params = [_SIZES] + param_names = ["size"] + + def setup(self, size): + numpy_random.seed(_SEED) + numpy_random.random_sample(size) + numpy_random.standard_normal(size) + numpy_random.normal(1.0, 2.0, size) + numpy_random.randint(0, 1000, size) + numpy_random.poisson(10.0, size) + + def time_random_sample(self, size): + numpy_random.random_sample(size) + + def time_standard_normal(self, size): + numpy_random.standard_normal(size) + + def time_normal(self, size): + numpy_random.normal(1.0, 2.0, size) + + def time_randint(self, size): + numpy_random.randint(0, 1000, size) + + def time_poisson(self, size): + numpy_random.poisson(10.0, size) + + +class PatchedNumpyRandom: + """numpy.random functions while patched by mkl_random. + + Fails instead of timing stock NumPy when the patch does not take effect. + """ + + params = [_SIZES] + param_names = ["size"] + + def setup(self, size): + mkl_random.patch_numpy_random() + served_by = np.random.standard_normal.__module__ or "" + if not (mkl_random.is_patched() and served_by.startswith("mkl_random")): + mkl_random.restore_numpy_random() + raise RuntimeError( + "[mkl-patch] numpy.random is not served by mkl_random " + f"after patch_numpy_random() (got {served_by!r})" + ) + np.random.seed(_SEED) + np.random.random_sample(size) + np.random.standard_normal(size) + np.random.randint(0, 1000, size) + + def teardown(self, size): + mkl_random.restore_numpy_random() + + def time_random_sample(self, size): + np.random.random_sample(size) + + def time_standard_normal(self, size): + np.random.standard_normal(size) + + def time_randint(self, size): + np.random.randint(0, 1000, size) diff --git a/benchmarks/benchmarks/bench_memory.py b/benchmarks/benchmarks/bench_memory.py new file mode 100644 index 0000000..8bff249 --- /dev/null +++ b/benchmarks/benchmarks/bench_memory.py @@ -0,0 +1,75 @@ +# Copyright (c) 2026, Intel Corporation +# +# 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 Intel Corporation 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 OWNER 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. + +"""Peak-memory benchmarks. + +Peak RSS includes setup, so setup allocates nothing large. +""" + +from ._utils import _make_state + +_N = 10_000_000 +_REPEATS = 10_000 + + +class PeakMemFill: + """Peak RSS of large fills, to catch new temporary buffers.""" + + def setup(self): + self.rs = _make_state() + + def peakmem_standard_normal(self): + self.rs.standard_normal(_N) + + def peakmem_randint_int64_small(self): + self.rs.randint(0, 100, _N, dtype="int64") + + def peakmem_randint_uint64_wide(self): + self.rs.randint(0, 3 * 2**32, _N, dtype="uint64") + + def peakmem_randint_uint8(self): + self.rs.randint(0, 2**8, _N, dtype="uint8") + + def peakmem_noncentral_chisquare(self): + self.rs.noncentral_chisquare(3.0, 2.0, _N) + + def peakmem_zipf(self): + self.rs.zipf(2.0, _N) + + +class PeakMemRepeat: + """Peak RSS of repeated small calls, to catch per-call leaks.""" + + def setup(self): + self.rs = _make_state() + self.state = self.rs.get_state() + + def peakmem_set_state(self): + for _ in range(_REPEATS): + self.rs.set_state(self.state) + + def peakmem_logseries(self): + for _ in range(_REPEATS): + self.rs.logseries(0.9, 1_000) diff --git a/benchmarks/benchmarks/bench_methods.py b/benchmarks/benchmarks/bench_methods.py new file mode 100644 index 0000000..a036436 --- /dev/null +++ b/benchmarks/benchmarks/bench_methods.py @@ -0,0 +1,116 @@ +# Copyright (c) 2026, Intel Corporation +# +# 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 Intel Corporation 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 OWNER 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. + +"""Benchmarks for the ``method`` keyword of mkl_random.MKLRandomState. + +Method names must be exact keys of mkl_random's alias tables: an unknown name +silently falls back to the default method. +""" + +import numpy as np + +from ._utils import _SEED, _SIZES, _make_state + +_GAUSSIAN_METHODS = ["ICDF", "BoxMuller", "BoxMuller2"] +# lognormal accepts no BoxMuller2 +_LOGNORMAL_METHODS = ["ICDF", "BoxMuller"] +_POISSON_METHODS = ["POISNORM", "PTPE"] +# MKL's PTPE poisson method may select a different internal algorithm +# depending on lam (around lam = 27); these values exercise both regimes +_POISSON_LAMS = [10.0, 100.0] + + +# --------------------------------------------------------------------------- +# Gaussian +# --------------------------------------------------------------------------- + + +class GaussianMethods: + """standard_normal / normal for each Gaussian method.""" + + params = [_GAUSSIAN_METHODS, _SIZES] + param_names = ["method", "size"] + + def setup(self, method, size): + self.rs = _make_state() + self.rs.standard_normal(size, method=method) + self.rs.normal(1.0, 2.0, size, method=method) + + def time_standard_normal(self, method, size): + self.rs.standard_normal(size, method=method) + + def time_normal(self, method, size): + self.rs.normal(1.0, 2.0, size, method=method) + + +class LognormalMethods: + """lognormal for each supported method.""" + + params = [_LOGNORMAL_METHODS, _SIZES] + param_names = ["method", "size"] + + def setup(self, method, size): + self.rs = _make_state() + self.rs.lognormal(0.0, 1.0, size, method=method) + + def time_lognormal(self, method, size): + self.rs.lognormal(0.0, 1.0, size, method=method) + + +class MultinormalCholesky: + """multinormal_cholesky (4-D) for each Gaussian method.""" + + params = [_GAUSSIAN_METHODS, _SIZES] + param_names = ["method", "size"] + + def setup(self, method, size): + self.rs = _make_state() + rng = np.random.default_rng(_SEED) + a = rng.standard_normal((4, 4)) + self.mean = rng.standard_normal(4) + self.ch = np.linalg.cholesky(a @ a.T + 4.0 * np.eye(4)) + self.rs.multinormal_cholesky(self.mean, self.ch, size, method=method) + + def time_multinormal_cholesky(self, method, size): + self.rs.multinormal_cholesky(self.mean, self.ch, size, method=method) + + +# --------------------------------------------------------------------------- +# Poisson +# --------------------------------------------------------------------------- + + +class PoissonMethods: + """poisson for each method, below and above the PTPE switch point.""" + + params = [_POISSON_METHODS, _POISSON_LAMS, _SIZES] + param_names = ["method", "lam", "size"] + + def setup(self, method, lam, size): + self.rs = _make_state() + self.rs.poisson(lam, size, method=method) + + def time_poisson(self, method, lam, size): + self.rs.poisson(lam, size, method=method) diff --git a/benchmarks/benchmarks/bench_multivariate.py b/benchmarks/benchmarks/bench_multivariate.py new file mode 100644 index 0000000..0df4b29 --- /dev/null +++ b/benchmarks/benchmarks/bench_multivariate.py @@ -0,0 +1,58 @@ +# Copyright (c) 2026, Intel Corporation +# +# 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 Intel Corporation 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 OWNER 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. + +"""Benchmarks for multivariate distributions; *size* is the sample count.""" + +import numpy as np + +from ._utils import _SEED, _SIZES, _make_state + + +class Multivariate: + """multinomial (5 categories), multivariate_normal (3-D), dirichlet (4-D)""" + + params = [_SIZES] + param_names = ["size"] + + def setup(self, size): + self.rs = _make_state() + rng = np.random.default_rng(_SEED) + a = rng.standard_normal((3, 3)) + self.mean = rng.standard_normal(3) + self.cov = a @ a.T + 3.0 * np.eye(3) + self.pvals = np.full(5, 0.2) + self.alpha = np.array([1.0, 2.0, 3.0, 4.0]) + self.rs.multinomial(20, self.pvals, size) + self.rs.multivariate_normal(self.mean, self.cov, size) + self.rs.dirichlet(self.alpha, size) + + def time_multinomial(self, size): + self.rs.multinomial(20, self.pvals, size) + + def time_multivariate_normal(self, size): + self.rs.multivariate_normal(self.mean, self.cov, size) + + def time_dirichlet(self, size): + self.rs.dirichlet(self.alpha, size) diff --git a/benchmarks/benchmarks/bench_permutations.py b/benchmarks/benchmarks/bench_permutations.py new file mode 100644 index 0000000..f43b6b0 --- /dev/null +++ b/benchmarks/benchmarks/bench_permutations.py @@ -0,0 +1,116 @@ +# Copyright (c) 2026, Intel Corporation +# +# 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 Intel Corporation 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 OWNER 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. + +"""Benchmarks for shuffle, permutation and choice.""" + +import numpy as np + +from ._utils import _SEED, _SIZES, _make_state + +_ROW = 8 # elements per row of the 2-D inputs + +# layout -> input builder. The layouts select the shuffle paths: +# 1d_int64 memcpy swaps, specialized for pointer-sized items +# 1d_int32 memcpy swaps, generic item size +# 2d_c memcpy swaps of contiguous rows +# 2d_f buffered swaps (rows are not contiguous) +# list untyped swaps of Python objects +_SHUFFLE_INPUTS = { + "1d_int64": lambda n: np.arange(n, dtype=np.int64), + "1d_int32": lambda n: np.arange(n, dtype=np.int32), + "2d_c": lambda n: np.arange(n, dtype=np.float64).reshape(-1, _ROW), + "2d_f": lambda n: np.asfortranarray( + np.arange(n, dtype=np.float64).reshape(-1, _ROW) + ), + "list": lambda n: list(range(n)), +} + +_CHOICE_POPULATION = 100_000 +_CHOICE_SAMPLES = 10_000 +# case -> (replace, weighted) +_CHOICE_CASES = { + "replace": (True, False), + "replace_p": (True, True), + "no_replace": (False, False), + "no_replace_p": (False, True), +} + + +class Shuffle: + """shuffle in place; *size* is the total number of elements.""" + + params = [list(_SHUFFLE_INPUTS), _SIZES] + param_names = ["layout", "size"] + + def setup(self, layout, size): + self.rs = _make_state() + self.x = _SHUFFLE_INPUTS[layout](size) + self.rs.shuffle(self.x) + + def time_shuffle(self, layout, size): + self.rs.shuffle(self.x) + + +class Permutation: + """permutation of a range, a 1-D array and the rows of a 2-D array.""" + + params = [["int", "1d", "2d"], _SIZES] + param_names = ["kind", "size"] + + def setup(self, kind, size): + self.rs = _make_state() + if kind == "int": + self.x = size + elif kind == "1d": + self.x = np.arange(size, dtype=np.int64) + else: + self.x = np.arange(size, dtype=np.float64).reshape(-1, _ROW) + self.rs.permutation(self.x) + + def time_permutation(self, kind, size): + self.rs.permutation(self.x) + + +class Choice: + """choice of 10k indices from 100k, with/without replacement and weights.""" + + params = [list(_CHOICE_CASES)] + param_names = ["case"] + + def setup(self, case): + self.rs = _make_state() + self.replace, weighted = _CHOICE_CASES[case] + self.p = None + if weighted: + w = np.random.default_rng(_SEED).random(_CHOICE_POPULATION) + self.p = w / w.sum() + self.rs.choice( + _CHOICE_POPULATION, _CHOICE_SAMPLES, replace=self.replace, p=self.p + ) + + def time_choice(self, case): + self.rs.choice( + _CHOICE_POPULATION, _CHOICE_SAMPLES, replace=self.replace, p=self.p + ) diff --git a/benchmarks/requirements.txt b/benchmarks/requirements.txt new file mode 100644 index 0000000..a4d92cc --- /dev/null +++ b/benchmarks/requirements.txt @@ -0,0 +1 @@ +psutil diff --git a/pyproject.toml b/pyproject.toml index 7a468b5..763dcfe 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -80,6 +80,7 @@ readme = { file = "README.md", content-type = "text/markdown" } requires-python = ">=3.10,<3.15" [project.optional-dependencies] +benchmark = ["asv>=0.6", "psutil"] test = ["pytest"] [project.urls]