Coverage for /usr/lib/python3/dist-packages/scipy/stats/_qmc.py: 14%
562 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Quasi-Monte Carlo engines and helpers."""
2from __future__ import annotations
4import copy
5import math
6import numbers
7import os
8import warnings
9from abc import ABC, abstractmethod
10from functools import partial
11from typing import (
12 Callable,
13 ClassVar,
14 Literal,
15 overload,
16 TYPE_CHECKING,
17)
19import numpy as np
21if TYPE_CHECKING:
22 import numpy.typing as npt
23 from scipy._lib._util import (
24 DecimalNumber, GeneratorType, IntNumber, SeedType
25 )
27import scipy.stats as stats
28from scipy._lib._util import rng_integers, _rng_spawn
29from scipy.spatial import distance, Voronoi
30from scipy.special import gammainc
31from ._sobol import (
32 _initialize_v, _cscramble, _fill_p_cumulative, _draw, _fast_forward,
33 _categorize, _MAXDIM
34)
35from ._qmc_cy import (
36 _cy_wrapper_centered_discrepancy,
37 _cy_wrapper_wrap_around_discrepancy,
38 _cy_wrapper_mixture_discrepancy,
39 _cy_wrapper_l2_star_discrepancy,
40 _cy_wrapper_update_discrepancy,
41 _cy_van_der_corput_scrambled,
42 _cy_van_der_corput,
43)
46__all__ = ['scale', 'discrepancy', 'update_discrepancy',
47 'QMCEngine', 'Sobol', 'Halton', 'LatinHypercube', 'PoissonDisk',
48 'MultinomialQMC', 'MultivariateNormalQMC']
51@overload
52def check_random_state(seed: IntNumber | None = ...) -> np.random.Generator:
53 ...
56@overload
57def check_random_state(seed: GeneratorType) -> GeneratorType:
58 ...
61# Based on scipy._lib._util.check_random_state
62def check_random_state(seed=None):
63 """Turn `seed` into a `numpy.random.Generator` instance.
65 Parameters
66 ----------
67 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional # noqa
68 If `seed` is an int or None, a new `numpy.random.Generator` is
69 created using ``np.random.default_rng(seed)``.
70 If `seed` is already a ``Generator`` or ``RandomState`` instance, then
71 the provided instance is used.
73 Returns
74 -------
75 seed : {`numpy.random.Generator`, `numpy.random.RandomState`}
76 Random number generator.
78 """
79 if seed is None or isinstance(seed, (numbers.Integral, np.integer)):
80 return np.random.default_rng(seed)
81 elif isinstance(seed, (np.random.RandomState, np.random.Generator)):
82 return seed
83 else:
84 raise ValueError(f'{seed!r} cannot be used to seed a'
85 ' numpy.random.Generator instance')
88def scale(
89 sample: npt.ArrayLike,
90 l_bounds: npt.ArrayLike,
91 u_bounds: npt.ArrayLike,
92 *,
93 reverse: bool = False
94) -> np.ndarray:
95 r"""Sample scaling from unit hypercube to different bounds.
97 To convert a sample from :math:`[0, 1)` to :math:`[a, b), b>a`,
98 with :math:`a` the lower bounds and :math:`b` the upper bounds.
99 The following transformation is used:
101 .. math::
103 (b - a) \cdot \text{sample} + a
105 Parameters
106 ----------
107 sample : array_like (n, d)
108 Sample to scale.
109 l_bounds, u_bounds : array_like (d,)
110 Lower and upper bounds (resp. :math:`a`, :math:`b`) of transformed
111 data. If `reverse` is True, range of the original data to transform
112 to the unit hypercube.
113 reverse : bool, optional
114 Reverse the transformation from different bounds to the unit hypercube.
115 Default is False.
117 Returns
118 -------
119 sample : array_like (n, d)
120 Scaled sample.
122 Examples
123 --------
124 Transform 3 samples in the unit hypercube to bounds:
126 >>> from scipy.stats import qmc
127 >>> l_bounds = [-2, 0]
128 >>> u_bounds = [6, 5]
129 >>> sample = [[0.5 , 0.75],
130 ... [0.5 , 0.5],
131 ... [0.75, 0.25]]
132 >>> sample_scaled = qmc.scale(sample, l_bounds, u_bounds)
133 >>> sample_scaled
134 array([[2. , 3.75],
135 [2. , 2.5 ],
136 [4. , 1.25]])
138 And convert back to the unit hypercube:
140 >>> sample_ = qmc.scale(sample_scaled, l_bounds, u_bounds, reverse=True)
141 >>> sample_
142 array([[0.5 , 0.75],
143 [0.5 , 0.5 ],
144 [0.75, 0.25]])
146 """
147 sample = np.asarray(sample)
149 # Checking bounds and sample
150 if not sample.ndim == 2:
151 raise ValueError('Sample is not a 2D array')
153 lower, upper = _validate_bounds(
154 l_bounds=l_bounds, u_bounds=u_bounds, d=sample.shape[1]
155 )
157 if not reverse:
158 # Checking that sample is within the hypercube
159 if (sample.max() > 1.) or (sample.min() < 0.):
160 raise ValueError('Sample is not in unit hypercube')
162 return sample * (upper - lower) + lower
163 else:
164 # Checking that sample is within the bounds
165 if not (np.all(sample >= lower) and np.all(sample <= upper)):
166 raise ValueError('Sample is out of bounds')
168 return (sample - lower) / (upper - lower)
171def discrepancy(
172 sample: npt.ArrayLike,
173 *,
174 iterative: bool = False,
175 method: Literal["CD", "WD", "MD", "L2-star"] = "CD",
176 workers: IntNumber = 1) -> float:
177 """Discrepancy of a given sample.
179 Parameters
180 ----------
181 sample : array_like (n, d)
182 The sample to compute the discrepancy from.
183 iterative : bool, optional
184 Must be False if not using it for updating the discrepancy.
185 Default is False. Refer to the notes for more details.
186 method : str, optional
187 Type of discrepancy, can be ``CD``, ``WD``, ``MD`` or ``L2-star``.
188 Refer to the notes for more details. Default is ``CD``.
189 workers : int, optional
190 Number of workers to use for parallel processing. If -1 is given all
191 CPU threads are used. Default is 1.
193 Returns
194 -------
195 discrepancy : float
196 Discrepancy.
198 Notes
199 -----
200 The discrepancy is a uniformity criterion used to assess the space filling
201 of a number of samples in a hypercube. A discrepancy quantifies the
202 distance between the continuous uniform distribution on a hypercube and the
203 discrete uniform distribution on :math:`n` distinct sample points.
205 The lower the value is, the better the coverage of the parameter space is.
207 For a collection of subsets of the hypercube, the discrepancy is the
208 difference between the fraction of sample points in one of those
209 subsets and the volume of that subset. There are different definitions of
210 discrepancy corresponding to different collections of subsets. Some
211 versions take a root mean square difference over subsets instead of
212 a maximum.
214 A measure of uniformity is reasonable if it satisfies the following
215 criteria [1]_:
217 1. It is invariant under permuting factors and/or runs.
218 2. It is invariant under rotation of the coordinates.
219 3. It can measure not only uniformity of the sample over the hypercube,
220 but also the projection uniformity of the sample over non-empty
221 subset of lower dimension hypercubes.
222 4. There is some reasonable geometric meaning.
223 5. It is easy to compute.
224 6. It satisfies the Koksma-Hlawka-like inequality.
225 7. It is consistent with other criteria in experimental design.
227 Four methods are available:
229 * ``CD``: Centered Discrepancy - subspace involves a corner of the
230 hypercube
231 * ``WD``: Wrap-around Discrepancy - subspace can wrap around bounds
232 * ``MD``: Mixture Discrepancy - mix between CD/WD covering more criteria
233 * ``L2-star``: L2-star discrepancy - like CD BUT variant to rotation
235 See [2]_ for precise definitions of each method.
237 Lastly, using ``iterative=True``, it is possible to compute the
238 discrepancy as if we had :math:`n+1` samples. This is useful if we want
239 to add a point to a sampling and check the candidate which would give the
240 lowest discrepancy. Then you could just update the discrepancy with
241 each candidate using `update_discrepancy`. This method is faster than
242 computing the discrepancy for a large number of candidates.
244 References
245 ----------
246 .. [1] Fang et al. "Design and modeling for computer experiments".
247 Computer Science and Data Analysis Series, 2006.
248 .. [2] Zhou Y.-D. et al. "Mixture discrepancy for quasi-random point sets."
249 Journal of Complexity, 29 (3-4) , pp. 283-301, 2013.
250 .. [3] T. T. Warnock. "Computational investigations of low discrepancy
251 point sets." Applications of Number Theory to Numerical
252 Analysis, Academic Press, pp. 319-343, 1972.
254 Examples
255 --------
256 Calculate the quality of the sample using the discrepancy:
258 >>> import numpy as np
259 >>> from scipy.stats import qmc
260 >>> space = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])
261 >>> l_bounds = [0.5, 0.5]
262 >>> u_bounds = [6.5, 6.5]
263 >>> space = qmc.scale(space, l_bounds, u_bounds, reverse=True)
264 >>> space
265 array([[0.08333333, 0.41666667],
266 [0.25 , 0.91666667],
267 [0.41666667, 0.25 ],
268 [0.58333333, 0.75 ],
269 [0.75 , 0.08333333],
270 [0.91666667, 0.58333333]])
271 >>> qmc.discrepancy(space)
272 0.008142039609053464
274 We can also compute iteratively the ``CD`` discrepancy by using
275 ``iterative=True``.
277 >>> disc_init = qmc.discrepancy(space[:-1], iterative=True)
278 >>> disc_init
279 0.04769081147119336
280 >>> qmc.update_discrepancy(space[-1], space[:-1], disc_init)
281 0.008142039609053513
283 """
284 sample = np.asarray(sample, dtype=np.float64, order="C")
286 # Checking that sample is within the hypercube and 2D
287 if not sample.ndim == 2:
288 raise ValueError("Sample is not a 2D array")
290 if (sample.max() > 1.) or (sample.min() < 0.):
291 raise ValueError("Sample is not in unit hypercube")
293 workers = _validate_workers(workers)
295 methods = {
296 "CD": _cy_wrapper_centered_discrepancy,
297 "WD": _cy_wrapper_wrap_around_discrepancy,
298 "MD": _cy_wrapper_mixture_discrepancy,
299 "L2-star": _cy_wrapper_l2_star_discrepancy,
300 }
302 if method in methods:
303 return methods[method](sample, iterative, workers=workers)
304 else:
305 raise ValueError(f"{method!r} is not a valid method. It must be one of"
306 f" {set(methods)!r}")
309def update_discrepancy(
310 x_new: npt.ArrayLike,
311 sample: npt.ArrayLike,
312 initial_disc: DecimalNumber) -> float:
313 """Update the centered discrepancy with a new sample.
315 Parameters
316 ----------
317 x_new : array_like (1, d)
318 The new sample to add in `sample`.
319 sample : array_like (n, d)
320 The initial sample.
321 initial_disc : float
322 Centered discrepancy of the `sample`.
324 Returns
325 -------
326 discrepancy : float
327 Centered discrepancy of the sample composed of `x_new` and `sample`.
329 Examples
330 --------
331 We can also compute iteratively the discrepancy by using
332 ``iterative=True``.
334 >>> import numpy as np
335 >>> from scipy.stats import qmc
336 >>> space = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])
337 >>> l_bounds = [0.5, 0.5]
338 >>> u_bounds = [6.5, 6.5]
339 >>> space = qmc.scale(space, l_bounds, u_bounds, reverse=True)
340 >>> disc_init = qmc.discrepancy(space[:-1], iterative=True)
341 >>> disc_init
342 0.04769081147119336
343 >>> qmc.update_discrepancy(space[-1], space[:-1], disc_init)
344 0.008142039609053513
346 """
347 sample = np.asarray(sample, dtype=np.float64, order="C")
348 x_new = np.asarray(x_new, dtype=np.float64, order="C")
350 # Checking that sample is within the hypercube and 2D
351 if not sample.ndim == 2:
352 raise ValueError('Sample is not a 2D array')
354 if (sample.max() > 1.) or (sample.min() < 0.):
355 raise ValueError('Sample is not in unit hypercube')
357 # Checking that x_new is within the hypercube and 1D
358 if not x_new.ndim == 1:
359 raise ValueError('x_new is not a 1D array')
361 if not (np.all(x_new >= 0) and np.all(x_new <= 1)):
362 raise ValueError('x_new is not in unit hypercube')
364 if x_new.shape[0] != sample.shape[1]:
365 raise ValueError("x_new and sample must be broadcastable")
367 return _cy_wrapper_update_discrepancy(x_new, sample, initial_disc)
370def _perturb_discrepancy(sample: np.ndarray, i1: int, i2: int, k: int,
371 disc: float):
372 """Centered discrepancy after an elementary perturbation of a LHS.
374 An elementary perturbation consists of an exchange of coordinates between
375 two points: ``sample[i1, k] <-> sample[i2, k]``. By construction,
376 this operation conserves the LHS properties.
378 Parameters
379 ----------
380 sample : array_like (n, d)
381 The sample (before permutation) to compute the discrepancy from.
382 i1 : int
383 The first line of the elementary permutation.
384 i2 : int
385 The second line of the elementary permutation.
386 k : int
387 The column of the elementary permutation.
388 disc : float
389 Centered discrepancy of the design before permutation.
391 Returns
392 -------
393 discrepancy : float
394 Centered discrepancy of the design after permutation.
396 References
397 ----------
398 .. [1] Jin et al. "An efficient algorithm for constructing optimal design
399 of computer experiments", Journal of Statistical Planning and
400 Inference, 2005.
402 """
403 n = sample.shape[0]
405 z_ij = sample - 0.5
407 # Eq (19)
408 c_i1j = (1. / n ** 2.
409 * np.prod(0.5 * (2. + abs(z_ij[i1, :])
410 + abs(z_ij) - abs(z_ij[i1, :] - z_ij)), axis=1))
411 c_i2j = (1. / n ** 2.
412 * np.prod(0.5 * (2. + abs(z_ij[i2, :])
413 + abs(z_ij) - abs(z_ij[i2, :] - z_ij)), axis=1))
415 # Eq (20)
416 c_i1i1 = (1. / n ** 2 * np.prod(1 + abs(z_ij[i1, :]))
417 - 2. / n * np.prod(1. + 0.5 * abs(z_ij[i1, :])
418 - 0.5 * z_ij[i1, :] ** 2))
419 c_i2i2 = (1. / n ** 2 * np.prod(1 + abs(z_ij[i2, :]))
420 - 2. / n * np.prod(1. + 0.5 * abs(z_ij[i2, :])
421 - 0.5 * z_ij[i2, :] ** 2))
423 # Eq (22), typo in the article in the denominator i2 -> i1
424 num = (2 + abs(z_ij[i2, k]) + abs(z_ij[:, k])
425 - abs(z_ij[i2, k] - z_ij[:, k]))
426 denum = (2 + abs(z_ij[i1, k]) + abs(z_ij[:, k])
427 - abs(z_ij[i1, k] - z_ij[:, k]))
428 gamma = num / denum
430 # Eq (23)
431 c_p_i1j = gamma * c_i1j
432 # Eq (24)
433 c_p_i2j = c_i2j / gamma
435 alpha = (1 + abs(z_ij[i2, k])) / (1 + abs(z_ij[i1, k]))
436 beta = (2 - abs(z_ij[i2, k])) / (2 - abs(z_ij[i1, k]))
438 g_i1 = np.prod(1. + abs(z_ij[i1, :]))
439 g_i2 = np.prod(1. + abs(z_ij[i2, :]))
440 h_i1 = np.prod(1. + 0.5 * abs(z_ij[i1, :]) - 0.5 * (z_ij[i1, :] ** 2))
441 h_i2 = np.prod(1. + 0.5 * abs(z_ij[i2, :]) - 0.5 * (z_ij[i2, :] ** 2))
443 # Eq (25), typo in the article g is missing
444 c_p_i1i1 = ((g_i1 * alpha) / (n ** 2) - 2. * alpha * beta * h_i1 / n)
445 # Eq (26), typo in the article n ** 2
446 c_p_i2i2 = ((g_i2 / ((n ** 2) * alpha)) - (2. * h_i2 / (n * alpha * beta)))
448 # Eq (26)
449 sum_ = c_p_i1j - c_i1j + c_p_i2j - c_i2j
451 mask = np.ones(n, dtype=bool)
452 mask[[i1, i2]] = False
453 sum_ = sum(sum_[mask])
455 disc_ep = (disc + c_p_i1i1 - c_i1i1 + c_p_i2i2 - c_i2i2 + 2 * sum_)
457 return disc_ep
460def primes_from_2_to(n: int) -> np.ndarray:
461 """Prime numbers from 2 to *n*.
463 Parameters
464 ----------
465 n : int
466 Sup bound with ``n >= 6``.
468 Returns
469 -------
470 primes : list(int)
471 Primes in ``2 <= p < n``.
473 Notes
474 -----
475 Taken from [1]_ by P.T. Roy, written consent given on 23.04.2021
476 by the original author, Bruno Astrolino, for free use in SciPy under
477 the 3-clause BSD.
479 References
480 ----------
481 .. [1] `StackOverflow <https://stackoverflow.com/questions/2068372>`_.
483 """
484 sieve = np.ones(n // 3 + (n % 6 == 2), dtype=bool)
485 for i in range(1, int(n ** 0.5) // 3 + 1):
486 k = 3 * i + 1 | 1
487 sieve[k * k // 3::2 * k] = False
488 sieve[k * (k - 2 * (i & 1) + 4) // 3::2 * k] = False
489 return np.r_[2, 3, ((3 * np.nonzero(sieve)[0][1:] + 1) | 1)]
492def n_primes(n: IntNumber) -> list[int]:
493 """List of the n-first prime numbers.
495 Parameters
496 ----------
497 n : int
498 Number of prime numbers wanted.
500 Returns
501 -------
502 primes : list(int)
503 List of primes.
505 """
506 primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
507 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127,
508 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193,
509 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269,
510 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
511 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431,
512 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503,
513 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
514 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673,
515 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761,
516 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857,
517 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947,
518 953, 967, 971, 977, 983, 991, 997][:n] # type: ignore[misc]
520 if len(primes) < n:
521 big_number = 2000
522 while 'Not enough primes':
523 primes = primes_from_2_to(big_number)[:n] # type: ignore
524 if len(primes) == n:
525 break
526 big_number += 1000
528 return primes
531def _van_der_corput_permutations(
532 base: IntNumber, *, random_state: SeedType = None
533) -> np.ndarray:
534 """Permutations for scrambling a Van der Corput sequence.
536 Parameters
537 ----------
538 base : int
539 Base of the sequence.
540 random_state : {None, int, `numpy.random.Generator`}, optional
541 If `seed` is an int or None, a new `numpy.random.Generator` is
542 created using ``np.random.default_rng(seed)``.
543 If `seed` is already a ``Generator`` instance, then the provided
544 instance is used.
546 Returns
547 -------
548 permutations : array_like
549 Permutation indices.
551 Notes
552 -----
553 In Algorithm 1 of Owen 2017, a permutation of `np.arange(base)` is
554 created for each positive integer `k` such that `1 - base**-k < 1`
555 using floating-point arithmetic. For double precision floats, the
556 condition `1 - base**-k < 1` can also be written as `base**-k >
557 2**-54`, which makes it more apparent how many permutations we need
558 to create.
559 """
560 rng = check_random_state(random_state)
561 count = math.ceil(54 / math.log2(base)) - 1
562 permutations = np.repeat(np.arange(base)[None], count, axis=0)
563 for perm in permutations:
564 rng.shuffle(perm)
566 return permutations
569def van_der_corput(
570 n: IntNumber,
571 base: IntNumber = 2,
572 *,
573 start_index: IntNumber = 0,
574 scramble: bool = False,
575 permutations: npt.ArrayLike | None = None,
576 seed: SeedType = None,
577 workers: IntNumber = 1) -> np.ndarray:
578 """Van der Corput sequence.
580 Pseudo-random number generator based on a b-adic expansion.
582 Scrambling uses permutations of the remainders (see [1]_). Multiple
583 permutations are applied to construct a point. The sequence of
584 permutations has to be the same for all points of the sequence.
586 Parameters
587 ----------
588 n : int
589 Number of element of the sequence.
590 base : int, optional
591 Base of the sequence. Default is 2.
592 start_index : int, optional
593 Index to start the sequence from. Default is 0.
594 scramble : bool, optional
595 If True, use Owen scrambling. Otherwise no scrambling is done.
596 Default is True.
597 permutations : array_like, optional
598 Permutations used for scrambling.
599 seed : {None, int, `numpy.random.Generator`}, optional
600 If `seed` is an int or None, a new `numpy.random.Generator` is
601 created using ``np.random.default_rng(seed)``.
602 If `seed` is already a ``Generator`` instance, then the provided
603 instance is used.
604 workers : int, optional
605 Number of workers to use for parallel processing. If -1 is
606 given all CPU threads are used. Default is 1.
608 Returns
609 -------
610 sequence : list (n,)
611 Sequence of Van der Corput.
613 References
614 ----------
615 .. [1] A. B. Owen. "A randomized Halton algorithm in R",
616 :arxiv:`1706.02808`, 2017.
618 """
619 if base < 2:
620 raise ValueError("'base' must be at least 2")
622 if scramble:
623 if permutations is None:
624 permutations = _van_der_corput_permutations(
625 base=base, random_state=seed
626 )
627 else:
628 permutations = np.asarray(permutations)
630 return _cy_van_der_corput_scrambled(n, base, start_index,
631 permutations, workers)
633 else:
634 return _cy_van_der_corput(n, base, start_index, workers)
637class QMCEngine(ABC):
638 """A generic Quasi-Monte Carlo sampler class meant for subclassing.
640 QMCEngine is a base class to construct a specific Quasi-Monte Carlo
641 sampler. It cannot be used directly as a sampler.
643 Parameters
644 ----------
645 d : int
646 Dimension of the parameter space.
647 optimization : {None, "random-cd", "lloyd"}, optional
648 Whether to use an optimization scheme to improve the quality after
649 sampling. Note that this is a post-processing step that does not
650 guarantee that all properties of the sample will be conserved.
651 Default is None.
653 * ``random-cd``: random permutations of coordinates to lower the
654 centered discrepancy. The best sample based on the centered
655 discrepancy is constantly updated. Centered discrepancy-based
656 sampling shows better space-filling robustness toward 2D and 3D
657 subprojections compared to using other discrepancy measures.
658 * ``lloyd``: Perturb samples using a modified Lloyd-Max algorithm.
659 The process converges to equally spaced samples.
661 .. versionadded:: 1.10.0
662 seed : {None, int, `numpy.random.Generator`}, optional
663 If `seed` is an int or None, a new `numpy.random.Generator` is
664 created using ``np.random.default_rng(seed)``.
665 If `seed` is already a ``Generator`` instance, then the provided
666 instance is used.
668 Notes
669 -----
670 By convention samples are distributed over the half-open interval
671 ``[0, 1)``. Instances of the class can access the attributes: ``d`` for
672 the dimension; and ``rng`` for the random number generator (used for the
673 ``seed``).
675 **Subclassing**
677 When subclassing `QMCEngine` to create a new sampler, ``__init__`` and
678 ``random`` must be redefined.
680 * ``__init__(d, seed=None)``: at least fix the dimension. If the sampler
681 does not take advantage of a ``seed`` (deterministic methods like
682 Halton), this parameter can be omitted.
683 * ``_random(n, *, workers=1)``: draw ``n`` from the engine. ``workers``
684 is used for parallelism. See `Halton` for example.
686 Optionally, two other methods can be overwritten by subclasses:
688 * ``reset``: Reset the engine to its original state.
689 * ``fast_forward``: If the sequence is deterministic (like Halton
690 sequence), then ``fast_forward(n)`` is skipping the ``n`` first draw.
692 Examples
693 --------
694 To create a random sampler based on ``np.random.random``, we would do the
695 following:
697 >>> from scipy.stats import qmc
698 >>> class RandomEngine(qmc.QMCEngine):
699 ... def __init__(self, d, seed=None):
700 ... super().__init__(d=d, seed=seed)
701 ...
702 ...
703 ... def _random(self, n=1, *, workers=1):
704 ... return self.rng.random((n, self.d))
705 ...
706 ...
707 ... def reset(self):
708 ... super().__init__(d=self.d, seed=self.rng_seed)
709 ... return self
710 ...
711 ...
712 ... def fast_forward(self, n):
713 ... self.random(n)
714 ... return self
716 After subclassing `QMCEngine` to define the sampling strategy we want to
717 use, we can create an instance to sample from.
719 >>> engine = RandomEngine(2)
720 >>> engine.random(5)
721 array([[0.22733602, 0.31675834], # random
722 [0.79736546, 0.67625467],
723 [0.39110955, 0.33281393],
724 [0.59830875, 0.18673419],
725 [0.67275604, 0.94180287]])
727 We can also reset the state of the generator and resample again.
729 >>> _ = engine.reset()
730 >>> engine.random(5)
731 array([[0.22733602, 0.31675834], # random
732 [0.79736546, 0.67625467],
733 [0.39110955, 0.33281393],
734 [0.59830875, 0.18673419],
735 [0.67275604, 0.94180287]])
737 """
739 @abstractmethod
740 def __init__(
741 self,
742 d: IntNumber,
743 *,
744 optimization: Literal["random-cd", "lloyd"] | None = None,
745 seed: SeedType = None
746 ) -> None:
747 if not np.issubdtype(type(d), np.integer) or d < 0:
748 raise ValueError('d must be a non-negative integer value')
750 self.d = d
752 if isinstance(seed, np.random.Generator):
753 # Spawn a Generator that we can own and reset.
754 self.rng = _rng_spawn(seed, 1)[0]
755 else:
756 # Create our instance of Generator, does not need spawning
757 # Also catch RandomState which cannot be spawned
758 self.rng = check_random_state(seed)
759 self.rng_seed = copy.deepcopy(self.rng)
761 self.num_generated = 0
763 config = {
764 # random-cd
765 "n_nochange": 100,
766 "n_iters": 10_000,
767 "rng": self.rng,
769 # lloyd
770 "tol": 1e-5,
771 "maxiter": 10,
772 "qhull_options": None,
773 }
774 self.optimization_method = _select_optimizer(optimization, config)
776 @abstractmethod
777 def _random(
778 self, n: IntNumber = 1, *, workers: IntNumber = 1
779 ) -> np.ndarray:
780 ...
782 def random(
783 self, n: IntNumber = 1, *, workers: IntNumber = 1
784 ) -> np.ndarray:
785 """Draw `n` in the half-open interval ``[0, 1)``.
787 Parameters
788 ----------
789 n : int, optional
790 Number of samples to generate in the parameter space.
791 Default is 1.
792 workers : int, optional
793 Only supported with `Halton`.
794 Number of workers to use for parallel processing. If -1 is
795 given all CPU threads are used. Default is 1. It becomes faster
796 than one worker for `n` greater than :math:`10^3`.
798 Returns
799 -------
800 sample : array_like (n, d)
801 QMC sample.
803 """
804 sample = self._random(n, workers=workers)
805 if self.optimization_method is not None:
806 sample = self.optimization_method(sample)
808 self.num_generated += n
809 return sample
811 def integers(
812 self,
813 l_bounds: npt.ArrayLike,
814 *,
815 u_bounds: npt.ArrayLike | None = None,
816 n: IntNumber = 1,
817 endpoint: bool = False,
818 workers: IntNumber = 1
819 ) -> np.ndarray:
820 r"""
821 Draw `n` integers from `l_bounds` (inclusive) to `u_bounds`
822 (exclusive), or if endpoint=True, `l_bounds` (inclusive) to
823 `u_bounds` (inclusive).
825 Parameters
826 ----------
827 l_bounds : int or array-like of ints
828 Lowest (signed) integers to be drawn (unless ``u_bounds=None``,
829 in which case this parameter is 0 and this value is used for
830 `u_bounds`).
831 u_bounds : int or array-like of ints, optional
832 If provided, one above the largest (signed) integer to be drawn
833 (see above for behavior if ``u_bounds=None``).
834 If array-like, must contain integer values.
835 n : int, optional
836 Number of samples to generate in the parameter space.
837 Default is 1.
838 endpoint : bool, optional
839 If true, sample from the interval ``[l_bounds, u_bounds]`` instead
840 of the default ``[l_bounds, u_bounds)``. Defaults is False.
841 workers : int, optional
842 Number of workers to use for parallel processing. If -1 is
843 given all CPU threads are used. Only supported when using `Halton`
844 Default is 1.
846 Returns
847 -------
848 sample : array_like (n, d)
849 QMC sample.
851 Notes
852 -----
853 It is safe to just use the same ``[0, 1)`` to integer mapping
854 with QMC that you would use with MC. You still get unbiasedness,
855 a strong law of large numbers, an asymptotically infinite variance
856 reduction and a finite sample variance bound.
858 To convert a sample from :math:`[0, 1)` to :math:`[a, b), b>a`,
859 with :math:`a` the lower bounds and :math:`b` the upper bounds,
860 the following transformation is used:
862 .. math::
864 \text{floor}((b - a) \cdot \text{sample} + a)
866 """
867 if u_bounds is None:
868 u_bounds = l_bounds
869 l_bounds = 0
871 u_bounds = np.atleast_1d(u_bounds)
872 l_bounds = np.atleast_1d(l_bounds)
874 if endpoint:
875 u_bounds = u_bounds + 1
877 if (not np.issubdtype(l_bounds.dtype, np.integer) or
878 not np.issubdtype(u_bounds.dtype, np.integer)):
879 message = ("'u_bounds' and 'l_bounds' must be integers or"
880 " array-like of integers")
881 raise ValueError(message)
883 if isinstance(self, Halton):
884 sample = self.random(n=n, workers=workers)
885 else:
886 sample = self.random(n=n)
888 sample = scale(sample, l_bounds=l_bounds, u_bounds=u_bounds)
889 sample = np.floor(sample).astype(np.int64)
891 return sample
893 def reset(self) -> QMCEngine:
894 """Reset the engine to base state.
896 Returns
897 -------
898 engine : QMCEngine
899 Engine reset to its base state.
901 """
902 seed = copy.deepcopy(self.rng_seed)
903 self.rng = check_random_state(seed)
904 self.num_generated = 0
905 return self
907 def fast_forward(self, n: IntNumber) -> QMCEngine:
908 """Fast-forward the sequence by `n` positions.
910 Parameters
911 ----------
912 n : int
913 Number of points to skip in the sequence.
915 Returns
916 -------
917 engine : QMCEngine
918 Engine reset to its base state.
920 """
921 self.random(n=n)
922 return self
925class Halton(QMCEngine):
926 """Halton sequence.
928 Pseudo-random number generator that generalize the Van der Corput sequence
929 for multiple dimensions. The Halton sequence uses the base-two Van der
930 Corput sequence for the first dimension, base-three for its second and
931 base-:math:`n` for its n-dimension.
933 Parameters
934 ----------
935 d : int
936 Dimension of the parameter space.
937 scramble : bool, optional
938 If True, use Owen scrambling. Otherwise no scrambling is done.
939 Default is True.
940 optimization : {None, "random-cd", "lloyd"}, optional
941 Whether to use an optimization scheme to improve the quality after
942 sampling. Note that this is a post-processing step that does not
943 guarantee that all properties of the sample will be conserved.
944 Default is None.
946 * ``random-cd``: random permutations of coordinates to lower the
947 centered discrepancy. The best sample based on the centered
948 discrepancy is constantly updated. Centered discrepancy-based
949 sampling shows better space-filling robustness toward 2D and 3D
950 subprojections compared to using other discrepancy measures.
951 * ``lloyd``: Perturb samples using a modified Lloyd-Max algorithm.
952 The process converges to equally spaced samples.
954 .. versionadded:: 1.10.0
955 seed : {None, int, `numpy.random.Generator`}, optional
956 If `seed` is an int or None, a new `numpy.random.Generator` is
957 created using ``np.random.default_rng(seed)``.
958 If `seed` is already a ``Generator`` instance, then the provided
959 instance is used.
961 Notes
962 -----
963 The Halton sequence has severe striping artifacts for even modestly
964 large dimensions. These can be ameliorated by scrambling. Scrambling
965 also supports replication-based error estimates and extends
966 applicabiltiy to unbounded integrands.
968 References
969 ----------
970 .. [1] Halton, "On the efficiency of certain quasi-random sequences of
971 points in evaluating multi-dimensional integrals", Numerische
972 Mathematik, 1960.
973 .. [2] A. B. Owen. "A randomized Halton algorithm in R",
974 :arxiv:`1706.02808`, 2017.
976 Examples
977 --------
978 Generate samples from a low discrepancy sequence of Halton.
980 >>> from scipy.stats import qmc
981 >>> sampler = qmc.Halton(d=2, scramble=False)
982 >>> sample = sampler.random(n=5)
983 >>> sample
984 array([[0. , 0. ],
985 [0.5 , 0.33333333],
986 [0.25 , 0.66666667],
987 [0.75 , 0.11111111],
988 [0.125 , 0.44444444]])
990 Compute the quality of the sample using the discrepancy criterion.
992 >>> qmc.discrepancy(sample)
993 0.088893711419753
995 If some wants to continue an existing design, extra points can be obtained
996 by calling again `random`. Alternatively, you can skip some points like:
998 >>> _ = sampler.fast_forward(5)
999 >>> sample_continued = sampler.random(n=5)
1000 >>> sample_continued
1001 array([[0.3125 , 0.37037037],
1002 [0.8125 , 0.7037037 ],
1003 [0.1875 , 0.14814815],
1004 [0.6875 , 0.48148148],
1005 [0.4375 , 0.81481481]])
1007 Finally, samples can be scaled to bounds.
1009 >>> l_bounds = [0, 2]
1010 >>> u_bounds = [10, 5]
1011 >>> qmc.scale(sample_continued, l_bounds, u_bounds)
1012 array([[3.125 , 3.11111111],
1013 [8.125 , 4.11111111],
1014 [1.875 , 2.44444444],
1015 [6.875 , 3.44444444],
1016 [4.375 , 4.44444444]])
1018 """
1020 def __init__(
1021 self, d: IntNumber, *, scramble: bool = True,
1022 optimization: Literal["random-cd", "lloyd"] | None = None,
1023 seed: SeedType = None
1024 ) -> None:
1025 # Used in `scipy.integrate.qmc_quad`
1026 self._init_quad = {'d': d, 'scramble': True,
1027 'optimization': optimization}
1028 super().__init__(d=d, optimization=optimization, seed=seed)
1029 self.seed = seed
1031 # important to have ``type(bdim) == int`` for performance reason
1032 self.base = [int(bdim) for bdim in n_primes(d)]
1033 self.scramble = scramble
1035 self._initialize_permutations()
1037 def _initialize_permutations(self) -> None:
1038 """Initialize permutations for all Van der Corput sequences.
1040 Permutations are only needed for scrambling.
1041 """
1042 self._permutations: list = [None] * len(self.base)
1043 if self.scramble:
1044 for i, bdim in enumerate(self.base):
1045 permutations = _van_der_corput_permutations(
1046 base=bdim, random_state=self.rng
1047 )
1049 self._permutations[i] = permutations
1051 def _random(
1052 self, n: IntNumber = 1, *, workers: IntNumber = 1
1053 ) -> np.ndarray:
1054 """Draw `n` in the half-open interval ``[0, 1)``.
1056 Parameters
1057 ----------
1058 n : int, optional
1059 Number of samples to generate in the parameter space. Default is 1.
1060 workers : int, optional
1061 Number of workers to use for parallel processing. If -1 is
1062 given all CPU threads are used. Default is 1. It becomes faster
1063 than one worker for `n` greater than :math:`10^3`.
1065 Returns
1066 -------
1067 sample : array_like (n, d)
1068 QMC sample.
1070 """
1071 workers = _validate_workers(workers)
1072 # Generate a sample using a Van der Corput sequence per dimension.
1073 sample = [van_der_corput(n, bdim, start_index=self.num_generated,
1074 scramble=self.scramble,
1075 permutations=self._permutations[i],
1076 workers=workers)
1077 for i, bdim in enumerate(self.base)]
1079 return np.array(sample).T.reshape(n, self.d)
1082class LatinHypercube(QMCEngine):
1083 r"""Latin hypercube sampling (LHS).
1085 A Latin hypercube sample [1]_ generates :math:`n` points in
1086 :math:`[0,1)^{d}`. Each univariate marginal distribution is stratified,
1087 placing exactly one point in :math:`[j/n, (j+1)/n)` for
1088 :math:`j=0,1,...,n-1`. They are still applicable when :math:`n << d`.
1090 Parameters
1091 ----------
1092 d : int
1093 Dimension of the parameter space.
1094 centered : bool, optional
1095 Center samples within cells of a multi-dimensional grid.
1096 Default is False.
1098 .. deprecated:: 1.10.0
1099 `centered` is deprecated as of SciPy 1.10.0 and will be removed in
1100 1.12.0. Use `scramble` instead. ``centered=True`` corresponds to
1101 ``scramble=False``.
1103 scramble : bool, optional
1104 When False, center samples within cells of a multi-dimensional grid.
1105 Otherwise, samples are randomly placed within cells of the grid.
1107 .. note::
1108 Setting ``scramble=False`` does not ensure deterministic output.
1109 For that, use the `seed` parameter.
1111 Default is True.
1113 .. versionadded:: 1.10.0
1115 optimization : {None, "random-cd", "lloyd"}, optional
1116 Whether to use an optimization scheme to improve the quality after
1117 sampling. Note that this is a post-processing step that does not
1118 guarantee that all properties of the sample will be conserved.
1119 Default is None.
1121 * ``random-cd``: random permutations of coordinates to lower the
1122 centered discrepancy. The best sample based on the centered
1123 discrepancy is constantly updated. Centered discrepancy-based
1124 sampling shows better space-filling robustness toward 2D and 3D
1125 subprojections compared to using other discrepancy measures.
1126 * ``lloyd``: Perturb samples using a modified Lloyd-Max algorithm.
1127 The process converges to equally spaced samples.
1129 .. versionadded:: 1.8.0
1130 .. versionchanged:: 1.10.0
1131 Add ``lloyd``.
1133 strength : {1, 2}, optional
1134 Strength of the LHS. ``strength=1`` produces a plain LHS while
1135 ``strength=2`` produces an orthogonal array based LHS of strength 2
1136 [7]_, [8]_. In that case, only ``n=p**2`` points can be sampled,
1137 with ``p`` a prime number. It also constrains ``d <= p + 1``.
1138 Default is 1.
1140 .. versionadded:: 1.8.0
1142 seed : {None, int, `numpy.random.Generator`}, optional
1143 If `seed` is an int or None, a new `numpy.random.Generator` is
1144 created using ``np.random.default_rng(seed)``.
1145 If `seed` is already a ``Generator`` instance, then the provided
1146 instance is used.
1148 Notes
1149 -----
1151 When LHS is used for integrating a function :math:`f` over :math:`n`,
1152 LHS is extremely effective on integrands that are nearly additive [2]_.
1153 With a LHS of :math:`n` points, the variance of the integral is always
1154 lower than plain MC on :math:`n-1` points [3]_. There is a central limit
1155 theorem for LHS on the mean and variance of the integral [4]_, but not
1156 necessarily for optimized LHS due to the randomization.
1158 :math:`A` is called an orthogonal array of strength :math:`t` if in each
1159 n-row-by-t-column submatrix of :math:`A`: all :math:`p^t` possible
1160 distinct rows occur the same number of times. The elements of :math:`A`
1161 are in the set :math:`\{0, 1, ..., p-1\}`, also called symbols.
1162 The constraint that :math:`p` must be a prime number is to allow modular
1163 arithmetic. Increasing strength adds some symmetry to the sub-projections
1164 of a sample. With strength 2, samples are symmetric along the diagonals of
1165 2D sub-projections. This may be undesirable, but on the other hand, the
1166 sample dispersion is improved.
1168 Strength 1 (plain LHS) brings an advantage over strength 0 (MC) and
1169 strength 2 is a useful increment over strength 1. Going to strength 3 is
1170 a smaller increment and scrambled QMC like Sobol', Halton are more
1171 performant [7]_.
1173 To create a LHS of strength 2, the orthogonal array :math:`A` is
1174 randomized by applying a random, bijective map of the set of symbols onto
1175 itself. For example, in column 0, all 0s might become 2; in column 1,
1176 all 0s might become 1, etc.
1177 Then, for each column :math:`i` and symbol :math:`j`, we add a plain,
1178 one-dimensional LHS of size :math:`p` to the subarray where
1179 :math:`A^i = j`. The resulting matrix is finally divided by :math:`p`.
1181 References
1182 ----------
1183 .. [1] Mckay et al., "A Comparison of Three Methods for Selecting Values
1184 of Input Variables in the Analysis of Output from a Computer Code."
1185 Technometrics, 1979.
1186 .. [2] M. Stein, "Large sample properties of simulations using Latin
1187 hypercube sampling." Technometrics 29, no. 2: 143-151, 1987.
1188 .. [3] A. B. Owen, "Monte Carlo variance of scrambled net quadrature."
1189 SIAM Journal on Numerical Analysis 34, no. 5: 1884-1910, 1997
1190 .. [4] Loh, W.-L. "On Latin hypercube sampling." The annals of statistics
1191 24, no. 5: 2058-2080, 1996.
1192 .. [5] Fang et al. "Design and modeling for computer experiments".
1193 Computer Science and Data Analysis Series, 2006.
1194 .. [6] Damblin et al., "Numerical studies of space filling designs:
1195 optimization of Latin Hypercube Samples and subprojection properties."
1196 Journal of Simulation, 2013.
1197 .. [7] A. B. Owen , "Orthogonal arrays for computer experiments,
1198 integration and visualization." Statistica Sinica, 1992.
1199 .. [8] B. Tang, "Orthogonal Array-Based Latin Hypercubes."
1200 Journal of the American Statistical Association, 1993.
1201 .. [9] Susan K. Seaholm et al. "Latin hypercube sampling and the
1202 sensitivity analysis of a Monte Carlo epidemic model".
1203 Int J Biomed Comput, 23(1-2), 97-112,
1204 :doi:`10.1016/0020-7101(88)90067-0`, 1988.
1206 Examples
1207 --------
1208 In [9]_, a Latin Hypercube sampling strategy was used to sample a
1209 parameter space to study the importance of each parameter of an epidemic
1210 model. Such analysis is also called a sensitivity analysis.
1212 Since the dimensionality of the problem is high (6), it is computationally
1213 expensive to cover the space. When numerical experiments are costly,
1214 QMC enables analysis that may not be possible if using a grid.
1216 The six parameters of the model represented the probability of illness,
1217 the probability of withdrawal, and four contact probabilities,
1218 The authors assumed uniform distributions for all parameters and generated
1219 50 samples.
1221 Using `scipy.stats.qmc.LatinHypercube` to replicate the protocol, the
1222 first step is to create a sample in the unit hypercube:
1224 >>> from scipy.stats import qmc
1225 >>> sampler = qmc.LatinHypercube(d=6)
1226 >>> sample = sampler.random(n=50)
1228 Then the sample can be scaled to the appropriate bounds:
1230 >>> l_bounds = [0.000125, 0.01, 0.0025, 0.05, 0.47, 0.7]
1231 >>> u_bounds = [0.000375, 0.03, 0.0075, 0.15, 0.87, 0.9]
1232 >>> sample_scaled = qmc.scale(sample, l_bounds, u_bounds)
1234 Such a sample was used to run the model 50 times, and a polynomial
1235 response surface was constructed. This allowed the authors to study the
1236 relative importance of each parameter across the range of
1237 possibilities of every other parameter.
1238 In this computer experiment, they showed a 14-fold reduction in the number
1239 of samples required to maintain an error below 2% on their response surface
1240 when compared to a grid sampling.
1242 Below are other examples showing alternative ways to construct LHS
1243 with even better coverage of the space.
1245 Using a base LHS as a baseline.
1247 >>> sampler = qmc.LatinHypercube(d=2)
1248 >>> sample = sampler.random(n=5)
1249 >>> qmc.discrepancy(sample)
1250 0.0196... # random
1252 Use the `optimization` keyword argument to produce a LHS with
1253 lower discrepancy at higher computational cost.
1255 >>> sampler = qmc.LatinHypercube(d=2, optimization="random-cd")
1256 >>> sample = sampler.random(n=5)
1257 >>> qmc.discrepancy(sample)
1258 0.0176... # random
1260 Use the `strength` keyword argument to produce an orthogonal array based
1261 LHS of strength 2. In this case, the number of sample points must be the
1262 square of a prime number.
1264 >>> sampler = qmc.LatinHypercube(d=2, strength=2)
1265 >>> sample = sampler.random(n=9)
1266 >>> qmc.discrepancy(sample)
1267 0.00526... # random
1269 Options could be combined to produce an optimized centered
1270 orthogonal array based LHS. After optimization, the result would not
1271 be guaranteed to be of strength 2.
1273 """
1275 def __init__(
1276 self, d: IntNumber, *, centered: bool = False,
1277 scramble: bool = True,
1278 strength: int = 1,
1279 optimization: Literal["random-cd", "lloyd"] | None = None,
1280 seed: SeedType = None
1281 ) -> None:
1282 if centered:
1283 scramble = False
1284 warnings.warn(
1285 "'centered' is deprecated and will be removed in SciPy 1.12."
1286 " Please use 'scramble' instead. 'centered=True' corresponds"
1287 " to 'scramble=False'.",
1288 stacklevel=2
1289 )
1291 # Used in `scipy.integrate.qmc_quad`
1292 self._init_quad = {'d': d, 'scramble': True, 'strength': strength,
1293 'optimization': optimization}
1294 super().__init__(d=d, seed=seed, optimization=optimization)
1295 self.scramble = scramble
1297 lhs_method_strength = {
1298 1: self._random_lhs,
1299 2: self._random_oa_lhs
1300 }
1302 try:
1303 self.lhs_method: Callable = lhs_method_strength[strength]
1304 except KeyError as exc:
1305 message = (f"{strength!r} is not a valid strength. It must be one"
1306 f" of {set(lhs_method_strength)!r}")
1307 raise ValueError(message) from exc
1309 def _random(
1310 self, n: IntNumber = 1, *, workers: IntNumber = 1
1311 ) -> np.ndarray:
1312 lhs = self.lhs_method(n)
1313 return lhs
1315 def _random_lhs(self, n: IntNumber = 1) -> np.ndarray:
1316 """Base LHS algorithm."""
1317 if not self.scramble:
1318 samples: np.ndarray | float = 0.5
1319 else:
1320 samples = self.rng.uniform(size=(n, self.d))
1322 perms = np.tile(np.arange(1, n + 1),
1323 (self.d, 1)) # type: ignore[arg-type]
1324 for i in range(self.d):
1325 self.rng.shuffle(perms[i, :])
1326 perms = perms.T
1328 samples = (perms - samples) / n
1329 return samples
1331 def _random_oa_lhs(self, n: IntNumber = 4) -> np.ndarray:
1332 """Orthogonal array based LHS of strength 2."""
1333 p = np.sqrt(n).astype(int)
1334 n_row = p**2
1335 n_col = p + 1
1337 primes = primes_from_2_to(p + 1)
1338 if p not in primes or n != n_row:
1339 raise ValueError(
1340 "n is not the square of a prime number. Close"
1341 f" values are {primes[-2:]**2}"
1342 )
1343 if self.d > p + 1:
1344 raise ValueError("n is too small for d. Must be n > (d-1)**2")
1346 oa_sample = np.zeros(shape=(n_row, n_col), dtype=int)
1348 # OA of strength 2
1349 arrays = np.tile(np.arange(p), (2, 1))
1350 oa_sample[:, :2] = np.stack(np.meshgrid(*arrays),
1351 axis=-1).reshape(-1, 2)
1352 for p_ in range(1, p):
1353 oa_sample[:, 2+p_-1] = np.mod(oa_sample[:, 0]
1354 + p_*oa_sample[:, 1], p)
1356 # scramble the OA
1357 oa_sample_ = np.empty(shape=(n_row, n_col), dtype=int)
1358 for j in range(n_col):
1359 perms = self.rng.permutation(p)
1360 oa_sample_[:, j] = perms[oa_sample[:, j]]
1362 # following is making a scrambled OA into an OA-LHS
1363 oa_lhs_sample = np.zeros(shape=(n_row, n_col))
1364 lhs_engine = LatinHypercube(d=1, scramble=self.scramble, strength=1,
1365 seed=self.rng) # type: QMCEngine
1366 for j in range(n_col):
1367 for k in range(p):
1368 idx = oa_sample[:, j] == k
1369 lhs = lhs_engine.random(p).flatten()
1370 oa_lhs_sample[:, j][idx] = lhs + oa_sample[:, j][idx]
1372 lhs_engine = lhs_engine.reset()
1374 oa_lhs_sample /= p
1376 return oa_lhs_sample[:, :self.d] # type: ignore
1379class Sobol(QMCEngine):
1380 """Engine for generating (scrambled) Sobol' sequences.
1382 Sobol' sequences are low-discrepancy, quasi-random numbers. Points
1383 can be drawn using two methods:
1385 * `random_base2`: safely draw :math:`n=2^m` points. This method
1386 guarantees the balance properties of the sequence.
1387 * `random`: draw an arbitrary number of points from the
1388 sequence. See warning below.
1390 Parameters
1391 ----------
1392 d : int
1393 Dimensionality of the sequence. Max dimensionality is 21201.
1394 scramble : bool, optional
1395 If True, use LMS+shift scrambling. Otherwise, no scrambling is done.
1396 Default is True.
1397 bits : int, optional
1398 Number of bits of the generator. Control the maximum number of points
1399 that can be generated, which is ``2**bits``. Maximal value is 64.
1400 It does not correspond to the return type, which is always
1401 ``np.float64`` to prevent points from repeating themselves.
1402 Default is None, which for backward compatibility, corresponds to 30.
1404 .. versionadded:: 1.9.0
1405 optimization : {None, "random-cd", "lloyd"}, optional
1406 Whether to use an optimization scheme to improve the quality after
1407 sampling. Note that this is a post-processing step that does not
1408 guarantee that all properties of the sample will be conserved.
1409 Default is None.
1411 * ``random-cd``: random permutations of coordinates to lower the
1412 centered discrepancy. The best sample based on the centered
1413 discrepancy is constantly updated. Centered discrepancy-based
1414 sampling shows better space-filling robustness toward 2D and 3D
1415 subprojections compared to using other discrepancy measures.
1416 * ``lloyd``: Perturb samples using a modified Lloyd-Max algorithm.
1417 The process converges to equally spaced samples.
1419 .. versionadded:: 1.10.0
1420 seed : {None, int, `numpy.random.Generator`}, optional
1421 If `seed` is an int or None, a new `numpy.random.Generator` is
1422 created using ``np.random.default_rng(seed)``.
1423 If `seed` is already a ``Generator`` instance, then the provided
1424 instance is used.
1426 Notes
1427 -----
1428 Sobol' sequences [1]_ provide :math:`n=2^m` low discrepancy points in
1429 :math:`[0,1)^{d}`. Scrambling them [3]_ makes them suitable for singular
1430 integrands, provides a means of error estimation, and can improve their
1431 rate of convergence. The scrambling strategy which is implemented is a
1432 (left) linear matrix scramble (LMS) followed by a digital random shift
1433 (LMS+shift) [2]_.
1435 There are many versions of Sobol' sequences depending on their
1436 'direction numbers'. This code uses direction numbers from [4]_. Hence,
1437 the maximum number of dimension is 21201. The direction numbers have been
1438 precomputed with search criterion 6 and can be retrieved at
1439 https://web.maths.unsw.edu.au/~fkuo/sobol/.
1441 .. warning::
1443 Sobol' sequences are a quadrature rule and they lose their balance
1444 properties if one uses a sample size that is not a power of 2, or skips
1445 the first point, or thins the sequence [5]_.
1447 If :math:`n=2^m` points are not enough then one should take :math:`2^M`
1448 points for :math:`M>m`. When scrambling, the number R of independent
1449 replicates does not have to be a power of 2.
1451 Sobol' sequences are generated to some number :math:`B` of bits.
1452 After :math:`2^B` points have been generated, the sequence would
1453 repeat. Hence, an error is raised.
1454 The number of bits can be controlled with the parameter `bits`.
1456 References
1457 ----------
1458 .. [1] I. M. Sobol', "The distribution of points in a cube and the accurate
1459 evaluation of integrals." Zh. Vychisl. Mat. i Mat. Phys., 7:784-802,
1460 1967.
1461 .. [2] J. Matousek, "On the L2-discrepancy for anchored boxes."
1462 J. of Complexity 14, 527-556, 1998.
1463 .. [3] Art B. Owen, "Scrambling Sobol and Niederreiter-Xing points."
1464 Journal of Complexity, 14(4):466-489, December 1998.
1465 .. [4] S. Joe and F. Y. Kuo, "Constructing sobol sequences with better
1466 two-dimensional projections." SIAM Journal on Scientific Computing,
1467 30(5):2635-2654, 2008.
1468 .. [5] Art B. Owen, "On dropping the first Sobol' point."
1469 :arxiv:`2008.08051`, 2020.
1471 Examples
1472 --------
1473 Generate samples from a low discrepancy sequence of Sobol'.
1475 >>> from scipy.stats import qmc
1476 >>> sampler = qmc.Sobol(d=2, scramble=False)
1477 >>> sample = sampler.random_base2(m=3)
1478 >>> sample
1479 array([[0. , 0. ],
1480 [0.5 , 0.5 ],
1481 [0.75 , 0.25 ],
1482 [0.25 , 0.75 ],
1483 [0.375, 0.375],
1484 [0.875, 0.875],
1485 [0.625, 0.125],
1486 [0.125, 0.625]])
1488 Compute the quality of the sample using the discrepancy criterion.
1490 >>> qmc.discrepancy(sample)
1491 0.013882107204860938
1493 To continue an existing design, extra points can be obtained
1494 by calling again `random_base2`. Alternatively, you can skip some
1495 points like:
1497 >>> _ = sampler.reset()
1498 >>> _ = sampler.fast_forward(4)
1499 >>> sample_continued = sampler.random_base2(m=2)
1500 >>> sample_continued
1501 array([[0.375, 0.375],
1502 [0.875, 0.875],
1503 [0.625, 0.125],
1504 [0.125, 0.625]])
1506 Finally, samples can be scaled to bounds.
1508 >>> l_bounds = [0, 2]
1509 >>> u_bounds = [10, 5]
1510 >>> qmc.scale(sample_continued, l_bounds, u_bounds)
1511 array([[3.75 , 3.125],
1512 [8.75 , 4.625],
1513 [6.25 , 2.375],
1514 [1.25 , 3.875]])
1516 """
1518 MAXDIM: ClassVar[int] = _MAXDIM
1520 def __init__(
1521 self, d: IntNumber, *, scramble: bool = True,
1522 bits: IntNumber | None = None, seed: SeedType = None,
1523 optimization: Literal["random-cd", "lloyd"] | None = None
1524 ) -> None:
1525 # Used in `scipy.integrate.qmc_quad`
1526 self._init_quad = {'d': d, 'scramble': True, 'bits': bits,
1527 'optimization': optimization}
1529 super().__init__(d=d, optimization=optimization, seed=seed)
1530 if d > self.MAXDIM:
1531 raise ValueError(
1532 f"Maximum supported dimensionality is {self.MAXDIM}."
1533 )
1535 self.bits = bits
1536 self.dtype_i: type
1538 if self.bits is None:
1539 self.bits = 30
1541 if self.bits <= 32:
1542 self.dtype_i = np.uint32
1543 elif 32 < self.bits <= 64:
1544 self.dtype_i = np.uint64
1545 else:
1546 raise ValueError("Maximum supported 'bits' is 64")
1548 self.maxn = 2**self.bits
1550 # v is d x maxbit matrix
1551 self._sv: np.ndarray = np.zeros((d, self.bits), dtype=self.dtype_i)
1552 _initialize_v(self._sv, dim=d, bits=self.bits)
1554 if not scramble:
1555 self._shift: np.ndarray = np.zeros(d, dtype=self.dtype_i)
1556 else:
1557 # scramble self._shift and self._sv
1558 self._scramble()
1560 self._quasi = self._shift.copy()
1562 # normalization constant with the largest possible number
1563 # calculate in Python to not overflow int with 2**64
1564 self._scale = 1.0 / 2 ** self.bits
1566 self._first_point = (self._quasi * self._scale).reshape(1, -1)
1567 # explicit casting to float64
1568 self._first_point = self._first_point.astype(np.float64)
1570 def _scramble(self) -> None:
1571 """Scramble the sequence using LMS+shift."""
1572 # Generate shift vector
1573 self._shift = np.dot(
1574 rng_integers(self.rng, 2, size=(self.d, self.bits),
1575 dtype=self.dtype_i),
1576 2 ** np.arange(self.bits, dtype=self.dtype_i),
1577 )
1578 # Generate lower triangular matrices (stacked across dimensions)
1579 ltm = np.tril(rng_integers(self.rng, 2,
1580 size=(self.d, self.bits, self.bits),
1581 dtype=self.dtype_i))
1582 _cscramble(
1583 dim=self.d, bits=self.bits, # type: ignore[arg-type]
1584 ltm=ltm, sv=self._sv
1585 )
1587 def _random(
1588 self, n: IntNumber = 1, *, workers: IntNumber = 1
1589 ) -> np.ndarray:
1590 """Draw next point(s) in the Sobol' sequence.
1592 Parameters
1593 ----------
1594 n : int, optional
1595 Number of samples to generate in the parameter space. Default is 1.
1597 Returns
1598 -------
1599 sample : array_like (n, d)
1600 Sobol' sample.
1602 """
1603 sample: np.ndarray = np.empty((n, self.d), dtype=np.float64)
1605 if n == 0:
1606 return sample
1608 total_n = self.num_generated + n
1609 if total_n > self.maxn:
1610 msg = (
1611 f"At most 2**{self.bits}={self.maxn} distinct points can be "
1612 f"generated. {self.num_generated} points have been previously "
1613 f"generated, then: n={self.num_generated}+{n}={total_n}. "
1614 )
1615 if self.bits != 64:
1616 msg += "Consider increasing `bits`."
1617 raise ValueError(msg)
1619 if self.num_generated == 0:
1620 # verify n is 2**n
1621 if not (n & (n - 1) == 0):
1622 warnings.warn("The balance properties of Sobol' points require"
1623 " n to be a power of 2.", stacklevel=2)
1625 if n == 1:
1626 sample = self._first_point
1627 else:
1628 _draw(
1629 n=n - 1, num_gen=self.num_generated, dim=self.d,
1630 scale=self._scale, sv=self._sv, quasi=self._quasi,
1631 sample=sample
1632 )
1633 sample = np.concatenate(
1634 [self._first_point, sample]
1635 )[:n] # type: ignore[misc]
1636 else:
1637 _draw(
1638 n=n, num_gen=self.num_generated - 1, dim=self.d,
1639 scale=self._scale, sv=self._sv, quasi=self._quasi,
1640 sample=sample
1641 )
1643 return sample
1645 def random_base2(self, m: IntNumber) -> np.ndarray:
1646 """Draw point(s) from the Sobol' sequence.
1648 This function draws :math:`n=2^m` points in the parameter space
1649 ensuring the balance properties of the sequence.
1651 Parameters
1652 ----------
1653 m : int
1654 Logarithm in base 2 of the number of samples; i.e., n = 2^m.
1656 Returns
1657 -------
1658 sample : array_like (n, d)
1659 Sobol' sample.
1661 """
1662 n = 2 ** m
1664 total_n = self.num_generated + n
1665 if not (total_n & (total_n - 1) == 0):
1666 raise ValueError("The balance properties of Sobol' points require "
1667 "n to be a power of 2. {0} points have been "
1668 "previously generated, then: n={0}+2**{1}={2}. "
1669 "If you still want to do this, the function "
1670 "'Sobol.random()' can be used."
1671 .format(self.num_generated, m, total_n))
1673 return self.random(n)
1675 def reset(self) -> Sobol:
1676 """Reset the engine to base state.
1678 Returns
1679 -------
1680 engine : Sobol
1681 Engine reset to its base state.
1683 """
1684 super().reset()
1685 self._quasi = self._shift.copy()
1686 return self
1688 def fast_forward(self, n: IntNumber) -> Sobol:
1689 """Fast-forward the sequence by `n` positions.
1691 Parameters
1692 ----------
1693 n : int
1694 Number of points to skip in the sequence.
1696 Returns
1697 -------
1698 engine : Sobol
1699 The fast-forwarded engine.
1701 """
1702 if self.num_generated == 0:
1703 _fast_forward(
1704 n=n - 1, num_gen=self.num_generated, dim=self.d,
1705 sv=self._sv, quasi=self._quasi
1706 )
1707 else:
1708 _fast_forward(
1709 n=n, num_gen=self.num_generated - 1, dim=self.d,
1710 sv=self._sv, quasi=self._quasi
1711 )
1712 self.num_generated += n
1713 return self
1716class PoissonDisk(QMCEngine):
1717 """Poisson disk sampling.
1719 Parameters
1720 ----------
1721 d : int
1722 Dimension of the parameter space.
1723 radius : float
1724 Minimal distance to keep between points when sampling new candidates.
1725 hypersphere : {"volume", "surface"}, optional
1726 Sampling strategy to generate potential candidates to be added in the
1727 final sample. Default is "volume".
1729 * ``volume``: original Bridson algorithm as described in [1]_.
1730 New candidates are sampled *within* the hypersphere.
1731 * ``surface``: only sample the surface of the hypersphere.
1732 ncandidates : int
1733 Number of candidates to sample per iteration. More candidates result
1734 in a denser sampling as more candidates can be accepted per iteration.
1735 optimization : {None, "random-cd", "lloyd"}, optional
1736 Whether to use an optimization scheme to improve the quality after
1737 sampling. Note that this is a post-processing step that does not
1738 guarantee that all properties of the sample will be conserved.
1739 Default is None.
1741 * ``random-cd``: random permutations of coordinates to lower the
1742 centered discrepancy. The best sample based on the centered
1743 discrepancy is constantly updated. Centered discrepancy-based
1744 sampling shows better space-filling robustness toward 2D and 3D
1745 subprojections compared to using other discrepancy measures.
1746 * ``lloyd``: Perturb samples using a modified Lloyd-Max algorithm.
1747 The process converges to equally spaced samples.
1749 .. versionadded:: 1.10.0
1750 seed : {None, int, `numpy.random.Generator`}, optional
1751 If `seed` is an int or None, a new `numpy.random.Generator` is
1752 created using ``np.random.default_rng(seed)``.
1753 If `seed` is already a ``Generator`` instance, then the provided
1754 instance is used.
1756 Notes
1757 -----
1758 Poisson disk sampling is an iterative sampling strategy. Starting from
1759 a seed sample, `ncandidates` are sampled in the hypersphere
1760 surrounding the seed. Candidates bellow a certain `radius` or outside the
1761 domain are rejected. New samples are added in a pool of sample seed. The
1762 process stops when the pool is empty or when the number of required
1763 samples is reached.
1765 The maximum number of point that a sample can contain is directly linked
1766 to the `radius`. As the dimension of the space increases, a higher radius
1767 spreads the points further and help overcome the curse of dimensionality.
1768 See the :ref:`quasi monte carlo tutorial <quasi-monte-carlo>` for more
1769 details.
1771 .. warning::
1773 The algorithm is more suitable for low dimensions and sampling size
1774 due to its iterative nature and memory requirements.
1775 Selecting a small radius with a high dimension would
1776 mean that the space could contain more samples than using lower
1777 dimension or a bigger radius.
1779 Some code taken from [2]_, written consent given on 31.03.2021
1780 by the original author, Shamis, for free use in SciPy under
1781 the 3-clause BSD.
1783 References
1784 ----------
1785 .. [1] Robert Bridson, "Fast Poisson Disk Sampling in Arbitrary
1786 Dimensions." SIGGRAPH, 2007.
1787 .. [2] `StackOverflow <https://stackoverflow.com/questions/66047540>`__.
1789 Examples
1790 --------
1791 Generate a 2D sample using a `radius` of 0.2.
1793 >>> import numpy as np
1794 >>> import matplotlib.pyplot as plt
1795 >>> from matplotlib.collections import PatchCollection
1796 >>> from scipy.stats import qmc
1797 >>>
1798 >>> rng = np.random.default_rng()
1799 >>> radius = 0.2
1800 >>> engine = qmc.PoissonDisk(d=2, radius=radius, seed=rng)
1801 >>> sample = engine.random(20)
1803 Visualizing the 2D sample and showing that no points are closer than
1804 `radius`. ``radius/2`` is used to visualize non-intersecting circles.
1805 If two samples are exactly at `radius` from each other, then their circle
1806 of radius ``radius/2`` will touch.
1808 >>> fig, ax = plt.subplots()
1809 >>> _ = ax.scatter(sample[:, 0], sample[:, 1])
1810 >>> circles = [plt.Circle((xi, yi), radius=radius/2, fill=False)
1811 ... for xi, yi in sample]
1812 >>> collection = PatchCollection(circles, match_original=True)
1813 >>> ax.add_collection(collection)
1814 >>> _ = ax.set(aspect='equal', xlabel=r'$x_1$', ylabel=r'$x_2$',
1815 ... xlim=[0, 1], ylim=[0, 1])
1816 >>> plt.show()
1818 Such visualization can be seen as circle packing: how many circle can
1819 we put in the space. It is a np-hard problem. The method `fill_space`
1820 can be used to add samples until no more samples can be added. This is
1821 a hard problem and parameters may need to be adjusted manually. Beware of
1822 the dimension: as the dimensionality increases, the number of samples
1823 required to fill the space increases exponentially
1824 (curse-of-dimensionality).
1826 """
1828 def __init__(
1829 self,
1830 d: IntNumber,
1831 *,
1832 radius: DecimalNumber = 0.05,
1833 hypersphere: Literal["volume", "surface"] = "volume",
1834 ncandidates: IntNumber = 30,
1835 optimization: Literal["random-cd", "lloyd"] | None = None,
1836 seed: SeedType = None
1837 ) -> None:
1838 # Used in `scipy.integrate.qmc_quad`
1839 self._init_quad = {'d': d, 'radius': radius,
1840 'hypersphere': hypersphere,
1841 'ncandidates': ncandidates,
1842 'optimization': optimization}
1843 super().__init__(d=d, optimization=optimization, seed=seed)
1845 hypersphere_sample = {
1846 "volume": self._hypersphere_volume_sample,
1847 "surface": self._hypersphere_surface_sample
1848 }
1850 try:
1851 self.hypersphere_method = hypersphere_sample[hypersphere]
1852 except KeyError as exc:
1853 message = (
1854 f"{hypersphere!r} is not a valid hypersphere sampling"
1855 f" method. It must be one of {set(hypersphere_sample)!r}")
1856 raise ValueError(message) from exc
1858 # size of the sphere from which the samples are drawn relative to the
1859 # size of a disk (radius)
1860 # for the surface sampler, all new points are almost exactly 1 radius
1861 # away from at least one existing sample +eps to avoid rejection
1862 self.radius_factor = 2 if hypersphere == "volume" else 1.001
1863 self.radius = radius
1864 self.radius_squared = self.radius**2
1866 # sample to generate per iteration in the hypersphere around center
1867 self.ncandidates = ncandidates
1869 with np.errstate(divide='ignore'):
1870 self.cell_size = self.radius / np.sqrt(self.d)
1871 self.grid_size = (
1872 np.ceil(np.ones(self.d) / self.cell_size)
1873 ).astype(int)
1875 self._initialize_grid_pool()
1877 def _initialize_grid_pool(self):
1878 """Sampling pool and sample grid."""
1879 self.sample_pool = []
1880 # Positions of cells
1881 # n-dim value for each grid cell
1882 self.sample_grid = np.empty(
1883 np.append(self.grid_size, self.d),
1884 dtype=np.float32
1885 )
1886 # Initialise empty cells with NaNs
1887 self.sample_grid.fill(np.nan)
1889 def _random(
1890 self, n: IntNumber = 1, *, workers: IntNumber = 1
1891 ) -> np.ndarray:
1892 """Draw `n` in the interval ``[0, 1]``.
1894 Note that it can return fewer samples if the space is full.
1895 See the note section of the class.
1897 Parameters
1898 ----------
1899 n : int, optional
1900 Number of samples to generate in the parameter space. Default is 1.
1902 Returns
1903 -------
1904 sample : array_like (n, d)
1905 QMC sample.
1907 """
1908 if n == 0 or self.d == 0:
1909 return np.empty((n, self.d))
1911 def in_limits(sample: np.ndarray) -> bool:
1912 return (sample.max() <= 1.) and (sample.min() >= 0.)
1914 def in_neighborhood(candidate: np.ndarray, n: int = 2) -> bool:
1915 """
1916 Check if there are samples closer than ``radius_squared`` to the
1917 `candidate` sample.
1918 """
1919 indices = (candidate / self.cell_size).astype(int)
1920 ind_min = np.maximum(indices - n, np.zeros(self.d, dtype=int))
1921 ind_max = np.minimum(indices + n + 1, self.grid_size)
1923 # Check if the center cell is empty
1924 if not np.isnan(self.sample_grid[tuple(indices)][0]):
1925 return True
1927 a = [slice(ind_min[i], ind_max[i]) for i in range(self.d)]
1929 # guards against: invalid value encountered in less as we are
1930 # comparing with nan and returns False. Which is wanted.
1931 with np.errstate(invalid='ignore'):
1932 if np.any(
1933 np.sum(
1934 np.square(candidate - self.sample_grid[tuple(a)]),
1935 axis=self.d
1936 ) < self.radius_squared
1937 ):
1938 return True
1940 return False
1942 def add_sample(candidate: np.ndarray) -> None:
1943 self.sample_pool.append(candidate)
1944 indices = (candidate / self.cell_size).astype(int)
1945 self.sample_grid[tuple(indices)] = candidate
1946 curr_sample.append(candidate)
1948 curr_sample: list[np.ndarray] = []
1950 if len(self.sample_pool) == 0:
1951 # the pool is being initialized with a single random sample
1952 add_sample(self.rng.random(self.d))
1953 num_drawn = 1
1954 else:
1955 num_drawn = 0
1957 # exhaust sample pool to have up to n sample
1958 while len(self.sample_pool) and num_drawn < n:
1959 # select a sample from the available pool
1960 idx_center = rng_integers(self.rng, len(self.sample_pool))
1961 center = self.sample_pool[idx_center]
1962 del self.sample_pool[idx_center]
1964 # generate candidates around the center sample
1965 candidates = self.hypersphere_method(
1966 center, self.radius * self.radius_factor, self.ncandidates
1967 )
1969 # keep candidates that satisfy some conditions
1970 for candidate in candidates:
1971 if in_limits(candidate) and not in_neighborhood(candidate):
1972 add_sample(candidate)
1974 num_drawn += 1
1975 if num_drawn >= n:
1976 break
1978 self.num_generated += num_drawn
1979 return np.array(curr_sample)
1981 def fill_space(self) -> np.ndarray:
1982 """Draw ``n`` samples in the interval ``[0, 1]``.
1984 Unlike `random`, this method will try to add points until
1985 the space is full. Depending on ``candidates`` (and to a lesser extent
1986 other parameters), some empty areas can still be present in the sample.
1988 .. warning::
1990 This can be extremely slow in high dimensions or if the
1991 ``radius`` is very small-with respect to the dimensionality.
1993 Returns
1994 -------
1995 sample : array_like (n, d)
1996 QMC sample.
1998 """
1999 return self.random(np.inf) # type: ignore[arg-type]
2001 def reset(self) -> PoissonDisk:
2002 """Reset the engine to base state.
2004 Returns
2005 -------
2006 engine : PoissonDisk
2007 Engine reset to its base state.
2009 """
2010 super().reset()
2011 self._initialize_grid_pool()
2012 return self
2014 def _hypersphere_volume_sample(
2015 self, center: np.ndarray, radius: DecimalNumber,
2016 candidates: IntNumber = 1
2017 ) -> np.ndarray:
2018 """Uniform sampling within hypersphere."""
2019 # should remove samples within r/2
2020 x = self.rng.standard_normal(size=(candidates, self.d))
2021 ssq = np.sum(x**2, axis=1)
2022 fr = radius * gammainc(self.d/2, ssq/2)**(1/self.d) / np.sqrt(ssq)
2023 fr_tiled = np.tile(
2024 fr.reshape(-1, 1), (1, self.d) # type: ignore[arg-type]
2025 )
2026 p = center + np.multiply(x, fr_tiled)
2027 return p
2029 def _hypersphere_surface_sample(
2030 self, center: np.ndarray, radius: DecimalNumber,
2031 candidates: IntNumber = 1
2032 ) -> np.ndarray:
2033 """Uniform sampling on the hypersphere's surface."""
2034 vec = self.rng.standard_normal(size=(candidates, self.d))
2035 vec /= np.linalg.norm(vec, axis=1)[:, None]
2036 p = center + np.multiply(vec, radius)
2037 return p
2040class MultivariateNormalQMC:
2041 r"""QMC sampling from a multivariate Normal :math:`N(\mu, \Sigma)`.
2043 Parameters
2044 ----------
2045 mean : array_like (d,)
2046 The mean vector. Where ``d`` is the dimension.
2047 cov : array_like (d, d), optional
2048 The covariance matrix. If omitted, use `cov_root` instead.
2049 If both `cov` and `cov_root` are omitted, use the identity matrix.
2050 cov_root : array_like (d, d'), optional
2051 A root decomposition of the covariance matrix, where ``d'`` may be less
2052 than ``d`` if the covariance is not full rank. If omitted, use `cov`.
2053 inv_transform : bool, optional
2054 If True, use inverse transform instead of Box-Muller. Default is True.
2055 engine : QMCEngine, optional
2056 Quasi-Monte Carlo engine sampler. If None, `Sobol` is used.
2057 seed : {None, int, `numpy.random.Generator`}, optional
2058 Used only if `engine` is None.
2059 If `seed` is an int or None, a new `numpy.random.Generator` is
2060 created using ``np.random.default_rng(seed)``.
2061 If `seed` is already a ``Generator`` instance, then the provided
2062 instance is used.
2064 Examples
2065 --------
2066 >>> import matplotlib.pyplot as plt
2067 >>> from scipy.stats import qmc
2068 >>> dist = qmc.MultivariateNormalQMC(mean=[0, 5], cov=[[1, 0], [0, 1]])
2069 >>> sample = dist.random(512)
2070 >>> _ = plt.scatter(sample[:, 0], sample[:, 1])
2071 >>> plt.show()
2073 """
2075 def __init__(
2076 self, mean: npt.ArrayLike, cov: npt.ArrayLike | None = None, *,
2077 cov_root: npt.ArrayLike | None = None,
2078 inv_transform: bool = True,
2079 engine: QMCEngine | None = None,
2080 seed: SeedType = None
2081 ) -> None:
2082 mean = np.array(mean, copy=False, ndmin=1)
2083 d = mean.shape[0]
2084 if cov is not None:
2085 # covariance matrix provided
2086 cov = np.array(cov, copy=False, ndmin=2)
2087 # check for square/symmetric cov matrix and mean vector has the
2088 # same d
2089 if not mean.shape[0] == cov.shape[0]:
2090 raise ValueError("Dimension mismatch between mean and "
2091 "covariance.")
2092 if not np.allclose(cov, cov.transpose()):
2093 raise ValueError("Covariance matrix is not symmetric.")
2094 # compute Cholesky decomp; if it fails, do the eigen decomposition
2095 try:
2096 cov_root = np.linalg.cholesky(cov).transpose()
2097 except np.linalg.LinAlgError:
2098 eigval, eigvec = np.linalg.eigh(cov)
2099 if not np.all(eigval >= -1.0e-8):
2100 raise ValueError("Covariance matrix not PSD.")
2101 eigval = np.clip(eigval, 0.0, None)
2102 cov_root = (eigvec * np.sqrt(eigval)).transpose()
2103 elif cov_root is not None:
2104 # root decomposition provided
2105 cov_root = np.atleast_2d(cov_root)
2106 if not mean.shape[0] == cov_root.shape[0]:
2107 raise ValueError("Dimension mismatch between mean and "
2108 "covariance.")
2109 else:
2110 # corresponds to identity covariance matrix
2111 cov_root = None
2113 self._inv_transform = inv_transform
2115 if not inv_transform:
2116 # to apply Box-Muller, we need an even number of dimensions
2117 engine_dim = 2 * math.ceil(d / 2)
2118 else:
2119 engine_dim = d
2120 if engine is None:
2121 self.engine = Sobol(
2122 d=engine_dim, scramble=True, bits=30, seed=seed
2123 ) # type: QMCEngine
2124 elif isinstance(engine, QMCEngine):
2125 if engine.d != engine_dim:
2126 raise ValueError("Dimension of `engine` must be consistent"
2127 " with dimensions of mean and covariance."
2128 " If `inv_transform` is False, it must be"
2129 " an even number.")
2130 self.engine = engine
2131 else:
2132 raise ValueError("`engine` must be an instance of "
2133 "`scipy.stats.qmc.QMCEngine` or `None`.")
2135 self._mean = mean
2136 self._corr_matrix = cov_root
2138 self._d = d
2140 def random(self, n: IntNumber = 1) -> np.ndarray:
2141 """Draw `n` QMC samples from the multivariate Normal.
2143 Parameters
2144 ----------
2145 n : int, optional
2146 Number of samples to generate in the parameter space. Default is 1.
2148 Returns
2149 -------
2150 sample : array_like (n, d)
2151 Sample.
2153 """
2154 base_samples = self._standard_normal_samples(n)
2155 return self._correlate(base_samples)
2157 def _correlate(self, base_samples: np.ndarray) -> np.ndarray:
2158 if self._corr_matrix is not None:
2159 return base_samples @ self._corr_matrix + self._mean
2160 else:
2161 # avoid multiplying with identity here
2162 return base_samples + self._mean
2164 def _standard_normal_samples(self, n: IntNumber = 1) -> np.ndarray:
2165 """Draw `n` QMC samples from the standard Normal :math:`N(0, I_d)`.
2167 Parameters
2168 ----------
2169 n : int, optional
2170 Number of samples to generate in the parameter space. Default is 1.
2172 Returns
2173 -------
2174 sample : array_like (n, d)
2175 Sample.
2177 """
2178 # get base samples
2179 samples = self.engine.random(n)
2180 if self._inv_transform:
2181 # apply inverse transform
2182 # (values to close to 0/1 result in inf values)
2183 return stats.norm.ppf(0.5 + (1 - 1e-10) * (samples - 0.5)) # type: ignore[attr-defined]
2184 else:
2185 # apply Box-Muller transform (note: indexes starting from 1)
2186 even = np.arange(0, samples.shape[-1], 2)
2187 Rs = np.sqrt(-2 * np.log(samples[:, even]))
2188 thetas = 2 * math.pi * samples[:, 1 + even]
2189 cos = np.cos(thetas)
2190 sin = np.sin(thetas)
2191 transf_samples = np.stack([Rs * cos, Rs * sin],
2192 -1).reshape(n, -1)
2193 # make sure we only return the number of dimension requested
2194 return transf_samples[:, : self._d]
2197class MultinomialQMC:
2198 r"""QMC sampling from a multinomial distribution.
2200 Parameters
2201 ----------
2202 pvals : array_like (k,)
2203 Vector of probabilities of size ``k``, where ``k`` is the number
2204 of categories. Elements must be non-negative and sum to 1.
2205 n_trials : int
2206 Number of trials.
2207 engine : QMCEngine, optional
2208 Quasi-Monte Carlo engine sampler. If None, `Sobol` is used.
2209 seed : {None, int, `numpy.random.Generator`}, optional
2210 Used only if `engine` is None.
2211 If `seed` is an int or None, a new `numpy.random.Generator` is
2212 created using ``np.random.default_rng(seed)``.
2213 If `seed` is already a ``Generator`` instance, then the provided
2214 instance is used.
2216 Examples
2217 --------
2218 Let's define 3 categories and for a given sample, the sum of the trials
2219 of each category is 8. The number of trials per category is determined
2220 by the `pvals` associated to each category.
2221 Then, we sample this distribution 64 times.
2223 >>> import matplotlib.pyplot as plt
2224 >>> from scipy.stats import qmc
2225 >>> dist = qmc.MultinomialQMC(
2226 ... pvals=[0.2, 0.4, 0.4], n_trials=10, engine=qmc.Halton(d=1)
2227 ... )
2228 >>> sample = dist.random(64)
2230 We can plot the sample and verify that the median of number of trials
2231 for each category is following the `pvals`. That would be
2232 ``pvals * n_trials = [2, 4, 4]``.
2234 >>> fig, ax = plt.subplots()
2235 >>> ax.yaxis.get_major_locator().set_params(integer=True)
2236 >>> _ = ax.boxplot(sample)
2237 >>> ax.set(xlabel="Categories", ylabel="Trials")
2238 >>> plt.show()
2240 """
2242 def __init__(
2243 self, pvals: npt.ArrayLike, n_trials: IntNumber,
2244 *, engine: QMCEngine | None = None,
2245 seed: SeedType = None
2246 ) -> None:
2247 self.pvals = np.array(pvals, copy=False, ndmin=1)
2248 if np.min(pvals) < 0:
2249 raise ValueError('Elements of pvals must be non-negative.')
2250 if not np.isclose(np.sum(pvals), 1):
2251 raise ValueError('Elements of pvals must sum to 1.')
2252 self.n_trials = n_trials
2253 if engine is None:
2254 self.engine = Sobol(
2255 d=1, scramble=True, bits=30, seed=seed
2256 ) # type: QMCEngine
2257 elif isinstance(engine, QMCEngine):
2258 if engine.d != 1:
2259 raise ValueError("Dimension of `engine` must be 1.")
2260 self.engine = engine
2261 else:
2262 raise ValueError("`engine` must be an instance of "
2263 "`scipy.stats.qmc.QMCEngine` or `None`.")
2265 def random(self, n: IntNumber = 1) -> np.ndarray:
2266 """Draw `n` QMC samples from the multinomial distribution.
2268 Parameters
2269 ----------
2270 n : int, optional
2271 Number of samples to generate in the parameter space. Default is 1.
2273 Returns
2274 -------
2275 samples : array_like (n, pvals)
2276 Sample.
2278 """
2279 sample = np.empty((n, len(self.pvals)))
2280 for i in range(n):
2281 base_draws = self.engine.random(self.n_trials).ravel()
2282 p_cumulative = np.empty_like(self.pvals, dtype=float)
2283 _fill_p_cumulative(np.array(self.pvals, dtype=float), p_cumulative)
2284 sample_ = np.zeros_like(self.pvals, dtype=int)
2285 _categorize(base_draws, p_cumulative, sample_)
2286 sample[i] = sample_
2287 return sample
2290def _select_optimizer(
2291 optimization: Literal["random-cd", "lloyd"] | None, config: dict
2292) -> Callable | None:
2293 """A factory for optimization methods."""
2294 optimization_method: dict[str, Callable] = {
2295 "random-cd": _random_cd,
2296 "lloyd": _lloyd_centroidal_voronoi_tessellation
2297 }
2299 optimizer: partial | None
2300 if optimization is not None:
2301 try:
2302 optimization = optimization.lower() # type: ignore[assignment]
2303 optimizer_ = optimization_method[optimization]
2304 except KeyError as exc:
2305 message = (f"{optimization!r} is not a valid optimization"
2306 f" method. It must be one of"
2307 f" {set(optimization_method)!r}")
2308 raise ValueError(message) from exc
2310 # config
2311 optimizer = partial(optimizer_, **config)
2312 else:
2313 optimizer = None
2315 return optimizer
2318def _random_cd(
2319 best_sample: np.ndarray, n_iters: int, n_nochange: int, rng: GeneratorType,
2320 **kwargs: dict
2321) -> np.ndarray:
2322 """Optimal LHS on CD.
2324 Create a base LHS and do random permutations of coordinates to
2325 lower the centered discrepancy.
2326 Because it starts with a normal LHS, it also works with the
2327 `centered` keyword argument.
2329 Two stopping criterion are used to stop the algorithm: at most,
2330 `n_iters` iterations are performed; or if there is no improvement
2331 for `n_nochange` consecutive iterations.
2332 """
2333 del kwargs # only use keywords which are defined, needed by factory
2335 n, d = best_sample.shape
2337 if d == 0 or n == 0:
2338 return np.empty((n, d))
2340 if d == 1 or n == 1:
2341 # discrepancy measures are invariant under permuting factors and runs
2342 return best_sample
2344 best_disc = discrepancy(best_sample)
2346 bounds = ([0, d - 1],
2347 [0, n - 1],
2348 [0, n - 1])
2350 n_nochange_ = 0
2351 n_iters_ = 0
2352 while n_nochange_ < n_nochange and n_iters_ < n_iters:
2353 n_iters_ += 1
2355 col = rng_integers(rng, *bounds[0], endpoint=True) # type: ignore[misc]
2356 row_1 = rng_integers(rng, *bounds[1], endpoint=True) # type: ignore[misc]
2357 row_2 = rng_integers(rng, *bounds[2], endpoint=True) # type: ignore[misc]
2358 disc = _perturb_discrepancy(best_sample,
2359 row_1, row_2, col,
2360 best_disc)
2361 if disc < best_disc:
2362 best_sample[row_1, col], best_sample[row_2, col] = (
2363 best_sample[row_2, col], best_sample[row_1, col])
2365 best_disc = disc
2366 n_nochange_ = 0
2367 else:
2368 n_nochange_ += 1
2370 return best_sample
2373def _l1_norm(sample: np.ndarray) -> float:
2374 return distance.pdist(sample, 'cityblock').min()
2377def _lloyd_iteration(
2378 sample: np.ndarray,
2379 decay: float,
2380 qhull_options: str
2381) -> np.ndarray:
2382 """Lloyd-Max algorithm iteration.
2384 Based on the implementation of Stéfan van der Walt:
2386 https://github.com/stefanv/lloyd
2388 which is:
2390 Copyright (c) 2021-04-21 Stéfan van der Walt
2391 https://github.com/stefanv/lloyd
2392 MIT License
2394 Parameters
2395 ----------
2396 sample : array_like (n, d)
2397 The sample to iterate on.
2398 decay : float
2399 Relaxation decay. A positive value would move the samples toward
2400 their centroid, and negative value would move them away.
2401 1 would move the samples to their centroid.
2402 qhull_options : str
2403 Additional options to pass to Qhull. See Qhull manual
2404 for details. (Default: "Qbb Qc Qz Qj Qx" for ndim > 4 and
2405 "Qbb Qc Qz Qj" otherwise.)
2407 Returns
2408 -------
2409 sample : array_like (n, d)
2410 The sample after an iteration of Lloyd's algorithm.
2412 """
2413 new_sample = np.empty_like(sample)
2415 voronoi = Voronoi(sample, qhull_options=qhull_options)
2417 for ii, idx in enumerate(voronoi.point_region):
2418 # the region is a series of indices into self.voronoi.vertices
2419 # remove samples at infinity, designated by index -1
2420 region = [i for i in voronoi.regions[idx] if i != -1]
2422 # get the vertices for this region
2423 verts = voronoi.vertices[region]
2425 # clipping would be wrong, we need to intersect
2426 # verts = np.clip(verts, 0, 1)
2428 # move samples towards centroids:
2429 # Centroid in n-D is the mean for uniformly distributed nodes
2430 # of a geometry.
2431 centroid = np.mean(verts, axis=0)
2432 new_sample[ii] = sample[ii] + (centroid - sample[ii]) * decay
2434 # only update sample to centroid within the region
2435 is_valid = np.all(np.logical_and(new_sample >= 0, new_sample <= 1), axis=1)
2436 sample[is_valid] = new_sample[is_valid]
2438 return sample
2441def _lloyd_centroidal_voronoi_tessellation(
2442 sample: npt.ArrayLike,
2443 *,
2444 tol: DecimalNumber = 1e-5,
2445 maxiter: IntNumber = 10,
2446 qhull_options: str | None = None,
2447 **kwargs: dict
2448) -> np.ndarray:
2449 """Approximate Centroidal Voronoi Tessellation.
2451 Perturb samples in N-dimensions using Lloyd-Max algorithm.
2453 Parameters
2454 ----------
2455 sample : array_like (n, d)
2456 The sample to iterate on. With ``n`` the number of samples and ``d``
2457 the dimension. Samples must be in :math:`[0, 1]^d`, with ``d>=2``.
2458 tol : float, optional
2459 Tolerance for termination. If the min of the L1-norm over the samples
2460 changes less than `tol`, it stops the algorithm. Default is 1e-5.
2461 maxiter : int, optional
2462 Maximum number of iterations. It will stop the algorithm even if
2463 `tol` is above the threshold.
2464 Too many iterations tend to cluster the samples as a hypersphere.
2465 Default is 10.
2466 qhull_options : str, optional
2467 Additional options to pass to Qhull. See Qhull manual
2468 for details. (Default: "Qbb Qc Qz Qj Qx" for ndim > 4 and
2469 "Qbb Qc Qz Qj" otherwise.)
2471 Returns
2472 -------
2473 sample : array_like (n, d)
2474 The sample after being processed by Lloyd-Max algorithm.
2476 Notes
2477 -----
2478 Lloyd-Max algorithm is an iterative process with the purpose of improving
2479 the dispersion of samples. For given sample: (i) compute a Voronoi
2480 Tessellation; (ii) find the centroid of each Voronoi cell; (iii) move the
2481 samples toward the centroid of their respective cell. See [1]_, [2]_.
2483 A relaxation factor is used to control how fast samples can move at each
2484 iteration. This factor is starting at 2 and ending at 1 after `maxiter`
2485 following an exponential decay.
2487 The process converges to equally spaced samples. It implies that measures
2488 like the discrepancy could suffer from too many iterations. On the other
2489 hand, L1 and L2 distances should improve. This is especially true with
2490 QMC methods which tend to favor the discrepancy over other criteria.
2492 .. note::
2494 The current implementation does not intersect the Voronoi Tessellation
2495 with the boundaries. This implies that for a low number of samples,
2496 empirically below 20, no Voronoi cell is touching the boundaries.
2497 Hence, samples cannot be moved close to the boundaries.
2499 Further improvements could consider the samples at infinity so that
2500 all boundaries are segments of some Voronoi cells. This would fix
2501 the computation of the centroid position.
2503 .. warning::
2505 The Voronoi Tessellation step is expensive and quickly becomes
2506 intractable with dimensions as low as 10 even for a sample
2507 of size as low as 1000.
2509 .. versionadded:: 1.9.0
2511 References
2512 ----------
2513 .. [1] Lloyd. "Least Squares Quantization in PCM".
2514 IEEE Transactions on Information Theory, 1982.
2515 .. [2] Max J. "Quantizing for minimum distortion".
2516 IEEE Transactions on Information Theory, 1960.
2518 Examples
2519 --------
2520 >>> import numpy as np
2521 >>> from scipy.spatial import distance
2522 >>> rng = np.random.default_rng()
2523 >>> sample = rng.random((128, 2))
2525 .. note::
2527 The samples need to be in :math:`[0, 1]^d`. `scipy.stats.qmc.scale`
2528 can be used to scale the samples from their
2529 original bounds to :math:`[0, 1]^d`. And back to their original bounds.
2531 Compute the quality of the sample using the L1 criterion.
2533 >>> def l1_norm(sample):
2534 ... return distance.pdist(sample, 'cityblock').min()
2536 >>> l1_norm(sample)
2537 0.00161... # random
2539 Now process the sample using Lloyd's algorithm and check the improvement
2540 on the L1. The value should increase.
2542 >>> sample = _lloyd_centroidal_voronoi_tessellation(sample)
2543 >>> l1_norm(sample)
2544 0.0278... # random
2546 """
2547 del kwargs # only use keywords which are defined, needed by factory
2549 sample = np.asarray(sample).copy()
2551 if not sample.ndim == 2:
2552 raise ValueError('`sample` is not a 2D array')
2554 if not sample.shape[1] >= 2:
2555 raise ValueError('`sample` dimension is not >= 2')
2557 # Checking that sample is within the hypercube
2558 if (sample.max() > 1.) or (sample.min() < 0.):
2559 raise ValueError('`sample` is not in unit hypercube')
2561 if qhull_options is None:
2562 qhull_options = 'Qbb Qc Qz QJ'
2564 if sample.shape[1] >= 5:
2565 qhull_options += ' Qx'
2567 # Fit an exponential to be 2 at 0 and 1 at `maxiter`.
2568 # The decay is used for relaxation.
2569 # analytical solution for y=exp(-maxiter/x) - 0.1
2570 root = -maxiter / np.log(0.1)
2571 decay = [np.exp(-x / root)+0.9 for x in range(maxiter)]
2573 l1_old = _l1_norm(sample=sample)
2574 for i in range(maxiter):
2575 sample = _lloyd_iteration(
2576 sample=sample, decay=decay[i],
2577 qhull_options=qhull_options,
2578 )
2580 l1_new = _l1_norm(sample=sample)
2582 if abs(l1_new - l1_old) < tol:
2583 break
2584 else:
2585 l1_old = l1_new
2587 return sample
2590def _validate_workers(workers: IntNumber = 1) -> IntNumber:
2591 """Validate `workers` based on platform and value.
2593 Parameters
2594 ----------
2595 workers : int, optional
2596 Number of workers to use for parallel processing. If -1 is
2597 given all CPU threads are used. Default is 1.
2599 Returns
2600 -------
2601 Workers : int
2602 Number of CPU used by the algorithm
2604 """
2605 workers = int(workers)
2606 if workers == -1:
2607 workers = os.cpu_count() # type: ignore[assignment]
2608 if workers is None:
2609 raise NotImplementedError(
2610 "Cannot determine the number of cpus using os.cpu_count(), "
2611 "cannot use -1 for the number of workers"
2612 )
2613 elif workers <= 0:
2614 raise ValueError(f"Invalid number of workers: {workers}, must be -1 "
2615 "or > 0")
2617 return workers
2620def _validate_bounds(
2621 l_bounds: npt.ArrayLike, u_bounds: npt.ArrayLike, d: int
2622) -> tuple[np.ndarray, ...]:
2623 """Bounds input validation.
2625 Parameters
2626 ----------
2627 l_bounds, u_bounds : array_like (d,)
2628 Lower and upper bounds.
2629 d : int
2630 Dimension to use for broadcasting.
2632 Returns
2633 -------
2634 l_bounds, u_bounds : array_like (d,)
2635 Lower and upper bounds.
2637 """
2638 try:
2639 lower = np.broadcast_to(l_bounds, d)
2640 upper = np.broadcast_to(u_bounds, d)
2641 except ValueError as exc:
2642 msg = ("'l_bounds' and 'u_bounds' must be broadcastable and respect"
2643 " the sample dimension")
2644 raise ValueError(msg) from exc
2646 if not np.all(lower < upper):
2647 raise ValueError("Bounds are not consistent 'l_bounds' < 'u_bounds'")
2649 return lower, upper