a a\ @s4dZddlmZddlZddlZddlZddlmZmZddl Z ddl m Z m Z m Z mZmZddlZddlZerddlm ZddlmZddlmZmZmZmZddlmZddlmZdd lm Z m!Z!m"Z"m#Z#m$Z$m%Z%m&Z&m'Z'm(Z(dd l)m*Z*m+Z+m,Z,m-Z-m.Z.gd Z/edCd ddddZ0edddddZ0dDddZ0ddddddddddZ1ddddddd d!d"d#d$d%Z2ddd&d"d'd(d)Z3d*dd+d,d-Z4d!d.d+d/d0Z5dEdddd2d!d!d!dd3dd4d5d6Z6Gd7d8d8eZ7Gd9d:d:e7Z8Gd;d<dd>e7Z:Gd?d@d@e7Z;GdAdBdBe7Zttjds2d}t|tj|St|tjj rP|St|tjj rb|Std|dS)aPTurn `seed` into a `numpy.random.Generator` instance. Parameters ---------- seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional If `seed` is None the `numpy.random.Generator` singleton is used. If `seed` is an int, a new ``Generator`` instance is used, seeded with `seed`. If `seed` is already a ``Generator`` or ``RandomState`` instance then that instance is used. Returns ------- seed : {`numpy.random.Generator`, `numpy.random.RandomState`} Random number generator. N GeneratorzaNumPy 1.16 doesn't have Generator, use either NumPy >= 1.17 or `seed=np.random.RandomState(seed)`z;%r cannot be used to seed a numpy.random.Generator instance) isinstancenumbersIntegralnpintegerhasattrrandom ValueErrorZ default_rngZ RandomStater.)r'msgr*r*r,r-:s  F)reverse npt.ArrayLikebool np.ndarray)samplel_boundsu_boundsr8r(cCst|}t|}t|}|jdks0tdt||\}}t||ksVtdt||jdkrptd|st|dkrt|dkstd||||St||krt||kstd||||Sd S) aSample scaling from unit hypercube to different bounds. To convert a sample from :math:`[0, 1)` to :math:`[a, b), b>a`, with :math:`a` the lower bounds and :math:`b` the upper bounds. The following transformation is used: .. math:: (b - a) \cdot \text{sample} + a Parameters ---------- sample : array_like (n, d) Sample to scale. l_bounds, u_bounds : array_like (d,) Lower and upper bounds (resp. :math:`a`, :math:`b`) of transformed data. If `reverse` is True, range of the original data to transform to the unit hypercube. reverse : bool, optional Reverse the transformation from different bounds to the unit hypercube. Default is False. Returns ------- sample : array_like (n, d) Scaled sample. Examples -------- Transform 3 samples in the unit hypercube to bounds: >>> from scipy.stats import qmc >>> l_bounds = [-2, 0] >>> u_bounds = [6, 5] >>> sample = [[0.5 , 0.75], ... [0.5 , 0.5], ... [0.75, 0.25]] >>> sample_scaled = qmc.scale(sample, l_bounds, u_bounds) >>> sample_scaled array([[2. , 3.75], [2. , 2.5 ], [4. , 1.25]]) And convert back to the unit hypercube: >>> sample_ = qmc.scale(sample_scaled, l_bounds, u_bounds, reverse=True) >>> sample_ array([[0.5 , 0.75], [0.5 , 0.5 ], [0.75, 0.25]]) Sample is not a 2D arrayzBounds are not consistent a < bz3Sample dimension is different than bounds dimensionrSample is not in unit hypercubezSample is out of boundsN) r2asarrayZ atleast_1dndimr6Zbroadcast_arraysalllenshape)r<r=r>r8lowerupperr*r*r,r_s";    rCDrA) iterativemethodworkersz$Literal['CD', 'WD', 'MD', 'L2-star']r float)r<rKrLrMr(cCstj|tjdd}|jdks$tdt|dkr@t|dksHtdt|}|dkrrt}|d urt d n|dkrtd |d t t t t d }||vr|||||dSt|dt|d S)aDiscrepancy of a given sample. Parameters ---------- sample : array_like (n, d) The sample to compute the discrepancy from. iterative : bool, optional Must be False if not using it for updating the discrepancy. Default is False. Refer to the notes for more details. method : str, optional Type of discrepancy, can be ``CD``, ``WD``, ``MD`` or ``L2-star``. Refer to the notes for more details. Default is ``CD``. workers : int, optional Number of workers to use for parallel processing. If -1 is given all CPU threads are used. Default is 1. Returns ------- discrepancy : float Discrepancy. Notes ----- The discrepancy is a uniformity criterion used to assess the space filling of a number of samples in a hypercube. A discrepancy quantifies the distance between the continuous uniform distribution on a hypercube and the discrete uniform distribution on :math:`n` distinct sample points. The lower the value is, the better the coverage of the parameter space is. For a collection of subsets of the hypercube, the discrepancy is the difference between the fraction of sample points in one of those subsets and the volume of that subset. There are different definitions of discrepancy corresponding to different collections of subsets. Some versions take a root mean square difference over subsets instead of a maximum. A measure of uniformity is reasonable if it satisfies the following criteria [1]_: 1. It is invariant under permuting factors and/or runs. 2. It is invariant under rotation of the coordinates. 3. It can measure not only uniformity of the sample over the hypercube, but also the projection uniformity of the sample over non-empty subset of lower dimension hypercubes. 4. There is some reasonable geometric meaning. 5. It is easy to compute. 6. It satisfies the Koksma-Hlawka-like inequality. 7. It is consistent with other criteria in experimental design. Four methods are available: * ``CD``: Centered Discrepancy - subspace involves a corner of the hypercube * ``WD``: Wrap-around Discrepancy - subspace can wrap around bounds * ``MD``: Mixture Discrepancy - mix between CD/WD covering more criteria * ``L2-star``: L2-star discrepancy - like CD BUT variant to rotation See [2]_ for precise definitions of each method. Lastly, using ``iterative=True``, it is possible to compute the discrepancy as if we had :math:`n+1` samples. This is useful if we want to add a point to a sampling and check the candidate which would give the lowest discrepancy. Then you could just update the discrepancy with each candidate using `update_discrepancy`. This method is faster than computing the discrepancy for a large number of candidates. References ---------- .. [1] Fang et al. "Design and modeling for computer experiments". Computer Science and Data Analysis Series, 2006. .. [2] Zhou Y.-D. et al. Mixture discrepancy for quasi-random point sets. Journal of Complexity, 29 (3-4) , pp. 283-301, 2013. .. [3] T. T. Warnock. "Computational investigations of low discrepancy point sets". Applications of Number Theory to Numerical Analysis, Academic Press, pp. 319-343, 1972. Examples -------- Calculate the quality of the sample using the discrepancy: >>> from scipy.stats import qmc >>> space = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]]) >>> l_bounds = [0.5, 0.5] >>> u_bounds = [6.5, 6.5] >>> space = qmc.scale(space, l_bounds, u_bounds, reverse=True) >>> space array([[0.08333333, 0.41666667], [0.25 , 0.91666667], [0.41666667, 0.25 ], [0.58333333, 0.75 ], [0.75 , 0.08333333], [0.91666667, 0.58333333]]) >>> qmc.discrepancy(space) 0.008142039609053464 We can also compute iteratively the ``CD`` discrepancy by using ``iterative=True``. >>> disc_init = qmc.discrepancy(space[:-1], iterative=True) >>> disc_init 0.04769081147119336 >>> qmc.update_discrepancy(space[-1], space[:-1], disc_init) 0.008142039609053513 Cdtypeorderr?r@rrArBNzaCannot determine the number of cpus using os.cpu_count(), cannot use -1 for the number of workerszInvalid number of workers: z, must be -1 or > 0)rJZWDZMDzL2-star)rMz* is not a valid method. It must be one of )r2rCfloat64rDr6rEintos cpu_countNotImplementedErrorrrrrset)r<rKrLrMmethodsr*r*r,rs0p rr )x_newr< initial_discr(cCstj|tjdd}tj|tjdd}|jdks6tdt|dkrRt|dksZtd|jdksltdt|dkrt|dkstd |jd|jdkrtd t|||S) aUpdate the centered discrepancy with a new sample. Parameters ---------- x_new : array_like (1, d) The new sample to add in `sample`. sample : array_like (n, d) The initial sample. initial_disc : float Centered discrepancy of the `sample`. Returns ------- discrepancy : float Centered discrepancy of the sample composed of `x_new` and `sample`. Examples -------- We can also compute iteratively the discrepancy by using ``iterative=True``. >>> from scipy.stats import qmc >>> space = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]]) >>> l_bounds = [0.5, 0.5] >>> u_bounds = [6.5, 6.5] >>> space = qmc.scale(space, l_bounds, u_bounds, reverse=True) >>> disc_init = qmc.discrepancy(space[:-1], iterative=True) >>> disc_init 0.04769081147119336 >>> qmc.update_discrepancy(space[-1], space[:-1], disc_init) 0.008142039609053513 rOrPr?r@rrArBzx_new is not a 1D arrayzx_new is not in unit hypercubez&x_new and sample must be broadcastable)r2rCrTrDr6rErGr)r[r<r\r*r*r,r Ks%  r rUnr(cCstj|d|ddktd}tdt|dddD]X}d|ddB}d|||ddd|<d|||d|d@d ddd|<q8tjdddt|d ddddBfS) aPrime numbers from 2 to *n*. Parameters ---------- n : int Sup bound with ``n >= 6``. Returns ------- primes : list(int) Primes in ``2 <= p < n``. Notes ----- Taken from [1]_ by P.T. Roy, written consent given on 23.04.2021 by the original author, Bruno Astrolino, for free use in SciPy under the 3-clause BSD. References ---------- .. [1] `StackOverflow `_. r?rQrA?FNr)r2Zonesr:rangerUZr_Znonzero)r^Zsieveikr*r*r,primes_from_2_tos ,rgz List[int]cCsLgdd|}t||krHd}t|d|}t||kr>qH|d7}q |S)zList of the n-first prime numbers. Parameters ---------- n : int Number of prime numbers wanted. Returns ------- primes : list(int) List of primes. )r?r_ %)+/5;=CGIOSYaegkmqiii iiiii%i3i7i9i=iKiQi[i]iaigioiui{iiiiiiiiiiiiiiiiiiiiiii i ii#i-i3i9i;iAiKiQiWiYi_ieiiikiwiiiiiiiiiiiiiiiiiiiiiiiii)i+i5i7i;i=iGiUiYi[i_imiqisiwiiiiiiiiiiiiiiNii)rFrg)r^ZprimesZ big_numberr*r*r,n_primess    rr? start_indexscrambler'r)r^baserrr'r(c Cst|}t|}t|||}d|}d|dkr||} |r`||} | t| t} || |7}||}|| |}q*|S)aVan der Corput sequence. Pseudo-random number generator based on a b-adic expansion. Scrambling uses permutations of the remainders (see [1]_). Multiple permutations are applied to construct a point. The sequence of permutations has to be the same for all points of the sequence. Parameters ---------- n : int Number of element of the sequence. base : int, optional Base of the sequence. Default is 2. start_index : int, optional Index to start the sequence from. Default is 0. scramble : bool, optional If True, use Owen scrambling. Otherwise no scrambling is done. Default is True. seed : {None, int, `numpy.random.Generator`}, optional If `seed` is None the `numpy.random.Generator` singleton is used. If `seed` is an int, a new ``Generator`` instance is used, seeded with `seed`. If `seed` is already a ``Generator`` instance then that instance is used. Returns ------- sequence : list (n,) Sequence of Van der Corput. References ---------- .. [1] A. B. Owen. "A randomized Halton algorithm in R", arXiv:1706.02808, 2017. rA)r-r2zerosarangeZ permutationarrayZastyperU) r^rrrr'rngsequenceZquotientZb2r remainderpermr*r*r,van_der_corputs,    rc@s`eZdZdZedddddddd Zeddd d d dZddddZddd ddZdS)r!a A generic Quasi-Monte Carlo sampler class meant for subclassing. QMCEngine is a base class to construct a specific Quasi-Monte Carlo sampler. It cannot be used directly as a sampler. Parameters ---------- d : int Dimension of the parameter space. seed : {None, int, `numpy.random.Generator`}, optional If `seed` is None the `numpy.random.Generator` singleton is used. If `seed` is an int, a new ``Generator`` instance is used, seeded with `seed`. If `seed` is already a ``Generator`` instance then that instance is used. Notes ----- By convention samples are distributed over the half-open interval ``[0, 1)``. Instances of the class can access the attributes: ``d`` for the dimension; and ``rng`` for the random number generator (used for the ``seed``). **Subclassing** When subclassing `QMCEngine` to create a new sampler, ``__init__`` and ``random`` must be redefined. * ``__init__(d, seed=None)``: at least fix the dimension. If the sampler does not take advantage of a ``seed`` (deterministic methods like Halton), this parameter can be omitted. * ``random(n)``: draw ``n`` from the engine and increase the counter ``num_generated`` by ``n``. Optionally, two other methods can be overwritten by subclasses: * ``reset``: Reset the engine to it's original state. * ``fast_forward``: If the sequence is deterministic (like Halton sequence), then ``fast_forward(n)`` is skipping the ``n`` first draw. Examples -------- To create a random sampler based on ``np.random.random``, we would do the following: >>> from scipy.stats import qmc >>> class RandomEngine(qmc.QMCEngine): ... def __init__(self, d, seed=None): ... super().__init__(d=d, seed=seed) ... ... ... def random(self, n=1): ... self.num_generated += n ... return self.rng.random((n, self.d)) ... ... ... def reset(self): ... super().__init__(d=self.d, seed=self.rng_seed) ... return self ... ... ... def fast_forward(self, n): ... self.random(n) ... return self After subclassing `QMCEngine` to define the sampling strategy we want to use, we can create an instance to sample from. >>> engine = RandomEngine(2) >>> engine.random(5) array([[0.22733602, 0.31675834], # random [0.79736546, 0.67625467], [0.39110955, 0.33281393], [0.59830875, 0.18673419], [0.67275604, 0.94180287]]) We can also reset the state of the generator and resample again. >>> _ = engine.reset() >>> engine.random(5) array([[0.22733602, 0.31675834], # random [0.79736546, 0.67625467], [0.39110955, 0.33281393], [0.59830875, 0.18673419], [0.67275604, 0.94180287]]) Nr+r rNone)dr'r(cCs@tt|tjstd||_t||_t ||_ d|_ dS)Nzd must be an integer valuer) r2Z issubdtypetyper3r6rr-rcopydeepcopyrng_seed num_generated)selfrr'r*r*r,__init__hs   zQMCEngine.__init__rAr;r]cCsdS)a1Draw `n` in the half-open interval ``[0, 1)``. Parameters ---------- n : int, optional Number of samples to generate in the parameter space. Default is 1. Returns ------- sample : array_like (n, d) QMC sample. Nr*rr^r*r*r,r5wszQMCEngine.randomr(cCs t|j}t||_d|_|S)zReset the engine to base state. Returns ------- engine : QMCEngine Engine reset to its base state. r)rrrr-rr)rr'r*r*r,resets  zQMCEngine.resetcCs|j|d|S)a Fast-forward the sequence by `n` positions. Parameters ---------- n : int Number of points to skip in the sequence. Returns ------- engine : QMCEngine Engine reset to its base state. )r^)r5rr*r*r, fast_forwards zQMCEngine.fast_forward)rA) __name__ __module__ __qualname____doc__rrr5rrr*r*r*r,r!sXr!csFeZdZdZdddddddd fd d Zddd dddZZS)r#a Halton sequence. Pseudo-random number generator that generalize the Van der Corput sequence for multiple dimensions. The Halton sequence uses the base-two Van der Corput sequence for the first dimension, base-three for its second and base-:math:`n` for its n-dimension. Parameters ---------- d : int Dimension of the parameter space. scramble : bool, optional If True, use Owen scrambling. Otherwise no scrambling is done. Default is True. seed : {None, int, `numpy.random.Generator`}, optional If `seed` is None the `numpy.random.Generator` singleton is used. If `seed` is an int, a new ``Generator`` instance is used, seeded with `seed`. If `seed` is already a ``Generator`` instance then that instance is used. Notes ----- The Halton sequence has severe striping artifacts for even modestly large dimensions. These can be ameliorated by scrambling. Scrambling also supports replication-based error estimates and extends applicabiltiy to unbounded integrands. References ---------- .. [1] Halton, "On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals", Numerische Mathematik, 1960. .. [2] A. B. Owen. "A randomized Halton algorithm in R", arXiv:1706.02808, 2017. Examples -------- Generate samples from a low discrepancy sequence of Halton. >>> from scipy.stats import qmc >>> sampler = qmc.Halton(d=2, scramble=False) >>> sample = sampler.random(n=5) >>> sample array([[0. , 0. ], [0.5 , 0.33333333], [0.25 , 0.66666667], [0.75 , 0.11111111], [0.125 , 0.44444444]]) Compute the quality of the sample using the discrepancy criterion. >>> qmc.discrepancy(sample) 0.088893711419753 If some wants to continue an existing design, extra points can be obtained by calling again `random`. Alternatively, you can skip some points like: >>> _ = sampler.fast_forward(5) >>> sample_continued = sampler.random(n=5) >>> sample_continued array([[0.3125 , 0.37037037], [0.8125 , 0.7037037 ], [0.1875 , 0.14814815], [0.6875 , 0.48148148], [0.4375 , 0.81481481]]) Finally, samples can be scaled to bounds. >>> l_bounds = [0, 2] >>> u_bounds = [10, 5] >>> qmc.scale(sample_continued, l_bounds, u_bounds) array([[3.125 , 3.11111111], [8.125 , 4.11111111], [1.875 , 2.44444444], [6.875 , 3.44444444], [4.375 , 4.44444444]]) TNrr'r r:rrrrr'r(cs*tj||d||_t||_||_dSNrr')superrr'rrrrrrr' __class__r*r,rs zHalton.__init__rAr;r]cs:fddjD}j7_t|jjS)a%Draw `n` in the half-open interval ``[0, 1)``. Parameters ---------- n : int, optional Number of samples to generate in the parameter space. Default is 1. Returns ------- sample : array_like (n, d) QMC sample. c s.g|]&}tt|jjtjdqS)r)rrUrrrrr').0Zbdimr^rr*r, s  z!Halton.random..)rrr2rTreshaperrr^r<r*rr,r5s  z Halton.random)rArrrrrr5 __classcell__r*r*rr,r#s Q r#csFeZdZdZdddddddd fd d Zddd dddZZS)r$a Latin hypercube sampling (LHS). A Latin hypercube sample [1]_ generates :math:`n` points in :math:`[0,1)^{d}`. Each univariate marginal distribution is stratified, placing exactly one point in :math:`[j/n, (j+1)/n)` for :math:`j=0,1,...,n-1`. They are still applicable when :math:`n << d`. LHS is extremely effective on integrands that are nearly additive [2]_. LHS on :math:`n` points never has more variance than plain MC on :math:`n-1` points [3]_. There is a central limit theorem for LHS [4]_, but not necessarily for optimized LHS. Parameters ---------- d : int Dimension of the parameter space. centered : bool, optional Center the point within the multi-dimensional grid. Default is False. seed : {None, int, `numpy.random.Generator`}, optional If `seed` is None the `numpy.random.Generator` singleton is used. If `seed` is an int, a new ``Generator`` instance is used, seeded with `seed`. If `seed` is already a ``Generator`` instance then that instance is used. References ---------- .. [1] Mckay et al., "A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code", Technometrics, 1979. .. [2] M. Stein, "Large sample properties of simulations using Latin hypercube sampling." Technometrics 29, no. 2: 143-151, 1987. .. [3] A. B. Owen, "Monte Carlo variance of scrambled net quadrature." SIAM Journal on Numerical Analysis 34, no. 5: 1884-1910, 1997 .. [4] Loh, W.-L. "On Latin hypercube sampling." The annals of statistics 24, no. 5: 2058-2080, 1996. Examples -------- Generate samples from a Latin hypercube generator. >>> from scipy.stats import qmc >>> sampler = qmc.LatinHypercube(d=2) >>> sample = sampler.random(n=5) >>> sample array([[0.1545328 , 0.53664833], # random [0.84052691, 0.06474907], [0.52177809, 0.93343721], [0.68033825, 0.36265316], [0.26544879, 0.61163943]]) Compute the quality of the sample using the discrepancy criterion. >>> qmc.discrepancy(sample) 0.019558034794794565 # random Finally, samples can be scaled to bounds. >>> l_bounds = [0, 2] >>> u_bounds = [10, 5] >>> qmc.scale(sample, l_bounds, u_bounds) array([[1.54532796, 3.609945 ], # random [8.40526909, 2.1942472 ], [5.2177809 , 4.80031164], [6.80338249, 3.08795949], [2.65448791, 3.83491828]]) FN)centeredr'r r:rr)rrr'r(cstj||d||_dSr)rrr)rrrr'rr*r,raszLatinHypercube.__init__rAr;r]cCs|jr d}n|jj||jfd}ttd|d|jdf}t|jD]}|j||ddfqH|j }|||}|j |7_ |S)a%Draw `n` in the half-open interval ``[0, 1)``. Parameters ---------- n : int, optional Number of samples to generate in the parameter space. Default is 1. Returns ------- sample : array_like (n, d) LHS sample. rb)sizerAN) rruniformrr2Ztilerrdshufflerr)rr^samplesZpermsrer*r*r,r5hs zLatinHypercube.random)rArr*r*rr,r$s Er$cseZdZUdZeZded<eZded<ddddd d d d fd dZ d dddZ ddddddZ dddddZ ddfdd Z dddddZZS)r"aEngine for generating (scrambled) Sobol' sequences. Sobol' sequences are low-discrepancy, quasi-random numbers. Points can be drawn using two methods: * `random_base2`: safely draw :math:`n=2^m` points. This method guarantees the balance properties of the sequence. * `random`: draw an arbitrary number of points from the sequence. See warning below. Parameters ---------- d : int Dimensionality of the sequence. Max dimensionality is 21201. scramble : bool, optional If True, use Owen scrambling. Otherwise no scrambling is done. Default is True. seed : {None, int, `numpy.random.Generator`}, optional If `seed` is None the `numpy.random.Generator` singleton is used. If `seed` is an int, a new ``Generator`` instance is used, seeded with `seed`. If `seed` is already a ``Generator`` instance then that instance is used. Notes ----- Sobol' sequences [1]_ provide :math:`n=2^m` low discrepancy points in :math:`[0,1)^{d}`. Scrambling them [2]_ makes them suitable for singular integrands, provides a means of error estimation, and can improve their rate of convergence. There are many versions of Sobol' sequences depending on their 'direction numbers'. This code uses direction numbers from [3]_. Hence, the maximum number of dimension is 21201. The direction numbers have been precomputed with search criterion 6 and can be retrieved at https://web.maths.unsw.edu.au/~fkuo/sobol/. .. warning:: Sobol' sequences are a quadrature rule and they lose their balance properties if one uses a sample size that is not a power of 2, or skips the first point, or thins the sequence [4]_. If :math:`n=2^m` points are not enough then one should take :math:`2^M` points for :math:`M>m`. When scrambling, the number R of independent replicates does not have to be a power of 2. Sobol' sequences are generated to some number :math:`B` of bits. After :math:`2^B` points have been generated, the sequence will repeat. Currently :math:`B=30`. References ---------- .. [1] I. M. Sobol. The distribution of points in a cube and the accurate evaluation of integrals. Zh. Vychisl. Mat. i Mat. Phys., 7:784-802, 1967. .. [2] Art B. Owen. Scrambling Sobol and Niederreiter-Xing points. Journal of Complexity, 14(4):466-489, December 1998. .. [3] S. Joe and F. Y. Kuo. Constructing sobol sequences with better two-dimensional projections. SIAM Journal on Scientific Computing, 30(5):2635-2654, 2008. .. [4] Art B. Owen. On dropping the first Sobol' point. arXiv 2008.08051, 2020. Examples -------- Generate samples from a low discrepancy sequence of Sobol'. >>> from scipy.stats import qmc >>> sampler = qmc.Sobol(d=2, scramble=False) >>> sample = sampler.random_base2(m=3) >>> sample array([[0. , 0. ], [0.5 , 0.5 ], [0.75 , 0.25 ], [0.25 , 0.75 ], [0.375, 0.375], [0.875, 0.875], [0.625, 0.125], [0.125, 0.625]]) Compute the quality of the sample using the discrepancy criterion. >>> qmc.discrepancy(sample) 0.013882107204860938 To continue an existing design, extra points can be obtained by calling again `random_base2`. Alternatively, you can skip some points like: >>> _ = sampler.reset() >>> _ = sampler.fast_forward(4) >>> sample_continued = sampler.random_base2(m=2) >>> sample_continued array([[0.375, 0.375], [0.875, 0.875], [0.625, 0.125], [0.125, 0.625]]) Finally, samples can be scaled to bounds. >>> l_bounds = [0, 2] >>> u_bounds = [10, 5] >>> qmc.scale(sample_continued, l_bounds, u_bounds) array([[3.75 , 3.125], [8.75 , 4.625], [6.25 , 2.375], [1.25 , 3.875]]) z ClassVar[int]MAXDIMMAXBITTNrr r:rrrcstj||d||jkr*td|jttj||jft d|_ t |j ||shtj|t d|_ n| |j |_|jd|jdd|_dS)Nrz'Maximum supported dimensionality is {}.rar?rArS)rrrr6formatrr2rrrU_svr_shift _scrambler_quasir _first_pointrrr*r,rs    zSobol.__init__rcCs~tt|jd|j|jftddtj|jtd|_|j |_ t t|jd|j|j|jftd}t |j||j d|_dS)zScramble the sequence.r?)rrQrarN)r2dotrrrrrUrrrrZtrilrrr)rZltmr*r*r,rs  zSobol._scramblerAr;r]cCstj||jftd}|jdkr||d@dks8td|dkrH|j}qt|d|j|j|j |j |t |j|gd|}nt||jd|j|j |j ||j|7_|S)a$Draw next point(s) in the Sobol' sequence. Parameters ---------- n : int, optional Number of samples to generate in the parameter space. Default is 1. Returns ------- sample : array_like (n, d) Sobol' sample. rarrAzEThe balance properties of Sobol' points require n to be a power of 2.N) r2emptyrrNrwarningswarnrrrrZ concatenaterr*r*r,r5#s  z Sobol.random)mr(cCs@d|}|j|}||d@dks6td|j||||S)aDraw point(s) from the Sobol' sequence. This function draws :math:`n=2^m` points in the parameter space ensuring the balance properties of the sequence. Parameters ---------- m : int Logarithm in base 2 of the number of samples; i.e., n = 2^m. Returns ------- sample : array_like (n, d) Sobol' sample. r?rArzThe balance properties of Sobol' points require n to be a power of 2. {0} points have been previously generated, then: n={0}+2**{1}={2}. If you still want to do this, the function 'Sobol.random()' can be used.)rr6rr5)rrr^Ztotal_nr*r*r, random_base2Fs zSobol.random_base2cst|j|_|S)zReset the engine to base state. Returns ------- engine : Sobol Engine reset to its base state. )rrrrrrrr*r,rds  z Sobol.resetcCsV|jdkr(t|d|j|j|j|jnt||jd|j|j|j|j|7_|S)aFast-forward the sequence by `n` positions. Parameters ---------- n : int Number of points to skip in the sequence. Returns ------- engine: Sobol The fast-forwarded engine. rrA)rrrrrrr*r*r,rqs zSobol.fast_forward)rA)rrrrrr__annotations__rrrrr5rrrrr*r*rr,r"s r  # r"c seZdZdZdddddddddddd d d fd d ZddddddZddfdd ZdddddZddddddZZ S)r&aiQMC sampling from a multivariate Normal :math:`N(\mu, \Sigma)`. Parameters ---------- mean : array_like (d,) The mean vector. Where ``d`` is the dimension. cov : array_like (d, d), optional The covariance matrix. If omitted, use `cov_root` instead. If both `cov` and `cov_root` are omitted, use the identity matrix. cov_root : array_like (d, d'), optional A root decomposition of the covariance matrix, where ``d'`` may be less than ``d`` if the covariance is not full rank. If omitted, use `cov`. inv_transform : bool, optional If True, use inverse transform instead of Box-Muller. Default is True. engine : QMCEngine, optional Quasi-Monte Carlo engine sampler. If None, `Sobol` is used. seed : {None, int, `numpy.random.Generator`}, optional If `seed` is None the `numpy.random.Generator` singleton is used. If `seed` is an int, a new ``Generator`` instance is used, seeded with `seed`. If `seed` is already a ``Generator`` instance then that instance is used. Examples -------- >>> import matplotlib.pyplot as plt >>> from scipy.stats import qmc >>> engine = qmc.MultivariateNormalQMC(mean=[0, 5], cov=[[1, 0], [0, 1]]) >>> sample = engine.random(512) >>> _ = plt.scatter(sample[:, 0], sample[:, 1]) >>> plt.show() NT)cov_root inv_transformenginer'r9zOptional[npt.ArrayLike]r:Optional[QMCEngine]rr)meancovrrrr'r(c stj|ddd}|jd}|durtj|ddd}|jd|jdksNtdt||sftdztj|}Wn\tjjytj |\}} t |dkstd t |d d}| t |}Yn0n8|durt |}|jd|jdkstdnd}tj||d ||_|sBdt|d} n|} |durbt| d |d |_n>> from scipy.stats import qmc >>> engine = qmc.MultinomialQMC(pvals=[0.2, 0.4, 0.4]) >>> sample = engine.random(10) N)rr'r9rrr)pvalsrr'r(cstj|ddd|_t|dkr(tdtt|dsBtd|dur\tdd|d|_n,t |t r|j dkrxtd ||_ntd t j d|d dS) NFrArrz'Elements of pvals must be non-negative.z Elements of pvals must sum to 1.Trz Dimension of `engine` must be 1.rr)r2rrminr6isclosesumr"rr/r!rrr)rrrr'rr*r,rHs  zMultinomialQMC.__init__rAr r;r]cCsd|j|}tj|jtd}ttj|jtd|tj |jt d}t ||||j |7_ |S)aJDraw `n` QMC samples from the multinomial distribution. Parameters ---------- n : int, optional Number of samples to generate in the parameter space. Default is 1. Returns ------- samples : array_like (pvals,) Vector of size ``p`` summing to `n`. ra) rr5Zravelr2Z empty_likerrNrrZ zeros_likerUrr)rr^Z base_drawsZ p_cumulativer<r*r*r,r5]s zMultinomialQMC.randomrcst|j|S)zReset the engine to base state. Returns ------- engine : MultinomialQMC Engine reset to its base state. rrrr*r,rss  zMultinomialQMC.reset)rA)rrrrrr5rrr*r*rr,r%/s r%).)N)r?)=r __future__rrVrr0abcrrrtypingrrrrr rZnumpyr2Z numpy.typingZnptZtyping_extensionsr Zscipy._lib._utilr r r rZ scipy.statsrrZscipy.stats._sobolrrrrrrrrrZscipy.stats._qmc_cyrrrrr__all__r-rrr rgrrr!r#r$r"r&r%r*r*r*r,sb     ,  *\< )Asi'