Discrete Distributions

UML Overview

`AbstractDiscreteDistribution`

Bases: object

Source code in qmcpy/discrete_distribution/abstract_discrete_distribution.py

def __init__(self, dimension, replications, seed, d_limit, n_limit):
    self.mimics = 'StdUniform'
    if not hasattr(self,'parameters'):
        self.parameters = []
    self.d_limit = d_limit
    self.n_limit = n_limit 
    if not (np.isscalar(self.d_limit) and self.d_limit>0 and np.isscalar(self.n_limit) and self.n_limit>0):
        raise ParameterError("d_limit and n_limit must be greater than 0")
    self.d_limit = np.inf if self.d_limit==np.inf else int(self.d_limit)
    self.n_limit = np.inf if self.n_limit==np.inf else int(self.n_limit)
    self.no_replications = replications is None 
    self.replications = 1 if self.no_replications else int(replications)
    if self.replications<0:
        raise ParameterError("replications must be None or a postive int")
    if isinstance(dimension,list) or isinstance(dimension,tuple) or isinstance(dimension,np.ndarray):
        self.dvec = np.array(dimension,dtype=int)
        self.d = len(self.dvec)
        if not (self.dvec.ndim==1 and len(np.unique(self.dvec))==self.d):
            raise ParameterError("dimension must be a 1d array of unique values")
    else:
        self.d = int(dimension)
        self.dvec = np.arange(self.d)
    if any(self.dvec>self.d_limit):
        raise ParameterError('dimension greater than dimension limit %d'%self.d_limit)
    self._base_seed = seed if isinstance(seed,np.random.SeedSequence) else np.random.SeedSequence(seed)
    self.entropy = self._base_seed.entropy
    self.spawn_key = self._base_seed.spawn_key
    self.rng = np.random.Generator(np.random.SFC64(self._base_seed))

call

__call__(n=None, n_min=None, n_max=None, return_binary=False, warn=True)

If just n is supplied, generate samples from the sequence at indices 0,...,n-1.
If n_min and n_max are supplied, generate samples from the sequence at indices n_min,...,n_max-1.
If n and n_min are supplied, then generate samples from the sequence at indices n,...,n_min-1.

Parameters:

Name	Type	Description	Default
`n`	`Union[None, int]`	Number of points to generate.	`None`
`n_min`	`Union[None, int]`	Starting index of sequence.	`None`
`n_max`	`Union[None, int]`	Final index of sequence.	`None`
`return_binary`	`bool`	Only used for `DigitalNetB2`. If `True`, only return the integer representation `x_integer` of base 2 digital net.	`False`
`warn`	`bool`	If `False`, disable warnings when generating samples.	`True`

Returns:

Name	Type	Description
`x`	`ndarray`	Samples from the sequence. If `replications` is `None` then this will be of size (`n_max`-`n_min`) \(\times\) `dimension` If `replications` is a positive int, then `x` will be of size `replications` \(\times\) (`n_max`-`n_min`) \(\times\) `dimension` Note that if `return_binary=True` then `x` is returned where `x` are integer representations of the digital net points.

Source code in qmcpy/discrete_distribution/abstract_discrete_distribution.py

def __call__(self, n=None, n_min=None, n_max=None, return_binary=False, warn=True):
    r"""
    - If just `n` is supplied, generate samples from the sequence at indices 0,...,`n`-1.
    - If `n_min` and `n_max` are supplied, generate samples from the sequence at indices `n_min`,...,`n_max`-1.
    - If `n` and `n_min` are supplied, then generate samples from the sequence at indices `n`,...,`n_min`-1.

    Args:
        n (Union[None,int]): Number of points to generate.
        n_min (Union[None,int]): Starting index of sequence.
        n_max (Union[None,int]): Final index of sequence.
        return_binary (bool): Only used for `DigitalNetB2`.  
            If `True`, *only* return the integer representation `x_integer` of base 2 digital net.  
        warn (bool): If `False`, disable warnings when generating samples.

    Returns:
        x (np.ndarray): Samples from the sequence. 

            - If `replications` is `None` then this will be of size (`n_max`-`n_min`) $\times$ `dimension` 
            - If `replications` is a positive int, then `x` will be of size `replications` $\times$ (`n_max`-`n_min`) $\times$ `dimension` 

            Note that if `return_binary=True` then `x` is returned where `x` are integer representations of the digital net points. 
    """
    return self.gen_samples(n=n,n_min=n_min,n_max=n_max,return_binary=return_binary,warn=warn)

spawn

spawn(s=1, dimensions=None)

Spawn new instances of the current discrete distribution but with new seeds and dimensions. Used by multi-level QMC algorithms which require different seeds and dimensions on each level.

Note

Use replications instead of using spawn when possible, e.g., when spawning copies which all have the same dimension.

Parameters:

Name	Type	Description	Default
`s`	`int`	Number of copies to spawn	`1`
`dimensions`	`ndarray`	Length `s` array of dimensions for each copy. Defaults to the current dimension.	`None`

Returns:

Name	Type	Description
`spawned_discrete_distribs`	`list`	Discrete distributions with new seeds and dimensions.

Source code in qmcpy/discrete_distribution/abstract_discrete_distribution.py

def spawn(self, s=1, dimensions=None):
    r"""
    Spawn new instances of the current discrete distribution but with new seeds and dimensions.
    Used by multi-level QMC algorithms which require different seeds and dimensions on each level.

    Note:
        Use `replications` instead of using `spawn` when possible, e.g., when spawning copies which all have the same dimension.

    Args:
        s (int): Number of copies to spawn
        dimensions (np.ndarray): Length `s` array of dimensions for each copy. Defaults to the current dimension. 

    Returns:
        spawned_discrete_distribs (list): Discrete distributions with new seeds and dimensions.
    """
    s = int(s)
    if s<=0:
        raise ParameterError("Must spawn s>0 instances")
    if dimensions is None: 
        dimensions = np.tile(self.d,s)
    elif (isinstance(dimensions,list) or isinstance(dimensions,tuple) or isinstance(dimensions,np.ndarray)):
        dimensions = np.array(dimensions,dtype=int)
    else:
        dimensions = np.tile(dimensions,s)
    if not (dimensions.ndim==1 and len(dimensions)==s):
        raise ParameterError("dimensions must be a length s np.ndarray")
    child_seeds = self._base_seed.spawn(s)
    spawned_discrete_distribs = [self._spawn(child_seeds[i],int(dimensions[i])) for i in range(s)]
    return spawned_discrete_distribs

`DigitalNetB2`

Bases: AbstractLDDiscreteDistribution

Low discrepancy digital net in base 2.

Note

Digital net sample sizes should be powers of \(2\) e.g. \(1\), \(2\), \(4\), \(8\), \(16\), \(\dots\).
The first point of an unrandomized digital nets is the origin.
Sobol is an alias for DigitalNetB2.
To use higher order digital nets, either:
- Pass in generating_matrices without interlacing and supply alpha>1 to apply interlacing, or
- Pass in generating_matrices with interlacing and set alpha=1 to avoid additional interlacing
i.e. do not pass in interlaced generating_matrices and set alpha>1, this will apply additional interlacing.

Examples:

>>> discrete_distrib = DigitalNetB2(2,seed=7)
>>> discrete_distrib(4)
array([[0.72162356, 0.914955  ],
       [0.16345554, 0.42964856],
       [0.98676255, 0.03436384],
       [0.42956655, 0.55876342]])
>>> discrete_distrib(1) # first point in the sequence
array([[0.72162356, 0.914955  ]])
>>> discrete_distrib
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               2^(1)
    replications    1
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         7

Replications of independent randomizations

>>> x = DigitalNetB2(dimension=3,seed=7,replications=2)(4)
>>> x.shape
(2, 4, 3)
>>> x
array([[[0.24653277, 0.1821862 , 0.74732591],
        [0.68152903, 0.66169442, 0.42891961],
        [0.48139855, 0.79818233, 0.08201287],
        [0.91541325, 0.29520621, 0.77495809]],

       [[0.44876891, 0.85899604, 0.50549679],
        [0.53635924, 0.04353443, 0.33564946],
        [0.23214143, 0.29281506, 0.06841036],
        [0.75295715, 0.60241448, 0.76962976]]])

Different orderings (avoid warnings that the first point is the origin)

>>> DigitalNetB2(dimension=2,randomize=False,order="GRAY")(n_min=2,n_max=4,warn=False)
array([[0.75, 0.25],
       [0.25, 0.75]])
>>> DigitalNetB2(dimension=2,randomize=False,order="RADICAL INVERSE")(n_min=2,n_max=4,warn=False)
array([[0.25, 0.75],
       [0.75, 0.25]])

Generating matrices from https://github.com/QMCSoftware/LDData/tree/main/dnet

>>> DigitalNetB2(dimension=3,randomize=False,generating_matrices="mps.nx_s5_alpha2_m32.txt")(8,warn=False)
array([[0.        , 0.        , 0.        ],
       [0.75841841, 0.45284834, 0.48844557],
       [0.57679828, 0.13226272, 0.10061957],
       [0.31858402, 0.32113875, 0.39369111],
       [0.90278927, 0.45867532, 0.01803333],
       [0.14542431, 0.02548793, 0.4749614 ],
       [0.45587539, 0.33081476, 0.11474426],
       [0.71318879, 0.15377192, 0.37629925]])

All randomizations

>>> DigitalNetB2(dimension=3,randomize='LMS DS',seed=5)(8)
array([[0.69346401, 0.20118185, 0.64779396],
       [0.43998032, 0.90102467, 0.0936172 ],
       [0.86663563, 0.60910036, 0.26043276],
       [0.11327376, 0.30772653, 0.93959283],
       [0.62102883, 0.79169756, 0.77051637],
       [0.37451038, 0.1231324 , 0.46634012],
       [0.94785596, 0.38577413, 0.13377215],
       [0.20121617, 0.71843325, 0.56293458]])
>>> DigitalNetB2(dimension=3,randomize='LMS',seed=5)(8,warn=False)
array([[0.        , 0.        , 0.        ],
       [0.75446077, 0.83265937, 0.69584079],
       [0.42329494, 0.65793842, 0.90427279],
       [0.67763292, 0.48937304, 0.33344964],
       [0.18550714, 0.97332905, 0.3772791 ],
       [0.93104851, 0.17195496, 0.82311652],
       [0.26221346, 0.31742386, 0.53093284],
       [0.50787715, 0.5172669 , 0.2101083 ]])
>>> DigitalNetB2(dimension=3,randomize='DS',seed=5)(8)
array([[0.68383949, 0.04047995, 0.42903182],
       [0.18383949, 0.54047995, 0.92903182],
       [0.93383949, 0.79047995, 0.67903182],
       [0.43383949, 0.29047995, 0.17903182],
       [0.55883949, 0.66547995, 0.05403182],
       [0.05883949, 0.16547995, 0.55403182],
       [0.80883949, 0.41547995, 0.80403182],
       [0.30883949, 0.91547995, 0.30403182]])
>>> DigitalNetB2(dimension=3,randomize='OWEN',seed=5)(8)
array([[0.33595486, 0.05834975, 0.30066401],
       [0.89110875, 0.84905188, 0.81833285],
       [0.06846074, 0.59997956, 0.67064205],
       [0.6693703 , 0.25824002, 0.10469644],
       [0.44586618, 0.99161977, 0.1873488 ],
       [0.84245267, 0.16445553, 0.56544372],
       [0.18546359, 0.44859876, 0.97389524],
       [0.61215442, 0.64341386, 0.44529863]])

Higher order net without randomization

>>> DigitalNetB2(dimension=3,randomize='FALSE',seed=7,alpha=2)(4,warn=False)
array([[0.    , 0.    , 0.    ],
       [0.75  , 0.75  , 0.75  ],
       [0.4375, 0.9375, 0.1875],
       [0.6875, 0.1875, 0.9375]])

Higher order nets with randomizations and replications

>>> DigitalNetB2(dimension=3,randomize='LMS DS',seed=7,replications=2,alpha=2)(4,warn=False)
array([[[0.42955149, 0.89149058, 0.43867111],
        [0.68701828, 0.07601148, 0.51312447],
        [0.10088033, 0.16293661, 0.25144138],
        [0.85846252, 0.87103178, 0.70041789]],

       [[0.27151905, 0.42406763, 0.21917369],
        [0.55035224, 0.67864387, 0.90033876],
        [0.19356758, 0.57589964, 0.00347701],
        [0.97235125, 0.32168581, 0.86920948]]])
>>> DigitalNetB2(dimension=3,randomize='LMS',seed=7,replications=2,alpha=2)(4,warn=False)
array([[[0.        , 0.        , 0.        ],
        [0.75817062, 0.96603053, 0.94947625],
        [0.45367986, 0.80295638, 0.18778553],
        [0.71171791, 0.2295424 , 0.76175441]],

       [[0.        , 0.        , 0.        ],
        [0.78664636, 0.75470215, 0.86876474],
        [0.45336727, 0.99953621, 0.22253579],
        [0.73996397, 0.24544824, 0.9008679 ]]])
>>> DigitalNetB2(dimension=3,randomize='DS',seed=7,replications=2,alpha=2)(4)
array([[[0.04386058, 0.58727432, 0.3691824 ],
        [0.79386058, 0.33727432, 0.6191824 ],
        [0.48136058, 0.39977432, 0.4316824 ],
        [0.73136058, 0.64977432, 0.6816824 ]],

       [[0.65212985, 0.69669968, 0.10605352],
        [0.40212985, 0.44669968, 0.85605352],
        [0.83962985, 0.25919968, 0.16855352],
        [0.08962985, 0.50919968, 0.91855352]]])
>>> DigitalNetB2(dimension=3,randomize='OWEN',seed=7,replications=2,alpha=2)(4)
array([[[0.46368517, 0.03964427, 0.62172094],
        [0.7498683 , 0.76141348, 0.4243043 ],
        [0.01729754, 0.97968459, 0.65963223],
        [0.75365329, 0.1903774 , 0.34141493]],

       [[0.52252547, 0.5679709 , 0.05949112],
        [0.27248656, 0.36488289, 0.81844058],
        [0.94219959, 0.39172304, 0.20285965],
        [0.19716391, 0.64741585, 0.92494554]]])

References:

Marius Hofert and Christiane Lemieux.
qrng: (Randomized) Quasi-Random Number Generators (2019).
R package version 0.0-7.
https://CRAN.R-project.org/package=qrng.
Faure, Henri, and Christiane Lemieux.
Implementation of Irreducible Sobol' Sequences in Prime Power Bases.
Mathematics and Computers in Simulation 161 (2019): 13-22. Crossref. Web.
F.Y. Kuo, D. Nuyens.
Application of quasi-Monte Carlo methods to elliptic PDEs with random diffusion coefficients - a survey of analysis and implementation.
Foundations of Computational Mathematics, 16(6):1631-1696, 2016.
https://link.springer.com/article/10.1007/s10208-016-9329-5.
D. Nuyens.
The Magic Point Shop of QMC point generators and generating vectors.
MATLAB and Python software, 2018.
https://people.cs.kuleuven.be/~dirk.nuyens/.
R. Cools, F.Y. Kuo, D. Nuyens.
Constructing embedded lattice rules for multivariate integration.
SIAM J. Sci. Comput., 28(6), 2162-2188.
I.M. Sobol', V.I. Turchaninov, Yu.L. Levitan, B.V. Shukhman.
Quasi-Random Sequence Generators.
Keldysh Institute of Applied Mathematics.
Russian Academy of Sciences, Moscow (1992).
Sobol, Ilya & Asotsky, Danil & Kreinin, Alexander & Kucherenko, Sergei. (2011).
Construction and Comparison of High-Dimensional Sobol' Generators. Wilmott. 2011.
10.1002/wilm.10056.
Paul Bratley and Bennett L. Fox.
Algorithm 659: Implementing Sobol's quasirandom sequence generator.
ACM Trans. Math. Softw. 14, 1 (March 1988), 88-100. 1988.
https://doi.org/10.1145/42288.214372.

Parameters:

Name	Type	Description	Default
`dimension`	`Union[int, ndarray]`	Dimension of the generator. If an `int` is passed in, use generating vector components at indices 0,...,`dimension`-1. If an `np.ndarray` is passed in, use generating vector components at these indices.	`1`
`replications`	`int`	Number of independent randomizations of a pointset.	`None`
`seed`	`Union[None,int,np.random.SeedSeq`	Seed the random number generator for reproducibility.	`None`
`randomize`	`str`	Options are `'LMS DS'`: Linear matrix scramble with digital shift. `'LMS'`: Linear matrix scramble only. `'DS'`: Digital shift only. `'NUS'`: Nested uniform scrambling. Also known as Owen scrambling. `'FALSE'`: No randomization. In this case the first point will be the origin.	`'LMS DS'`
`generating_matrices`	`Union[str, ndarray, int]`	Specify the generating matrices. A `str` should be the name (or path) of a file from the LDData repo at https://github.com/QMCSoftware/LDData/tree/main/dnet. An `np.ndarray` of integers with shape \((d,m_\mathrm{max})\) or \((r,d,m_\mathrm{max})\) where \(d\) is the number of dimensions, \(r\) is the number of replications, and \(2^{m_\mathrm{max}}\) is the maximum number of supported points. Setting `msb=False` will flip the bits of ints in the generating matrices.	`'joe_kuo.6.21201.txt'`
`order`	`str`	`'RADICAL INVERSE'`, or `'GRAY'` ordering. See the doctest example above.	`'RADICAL INVERSE'`
`t`	`int`	Number of bits in integer represetation of points after randomization. The number of bits in the generating matrices is inferred based on the largest value.	`63`
`alpha`	`int`	Interlacing factor for higher order nets. When `alpha`>1, interlacing is performed regardless of the generating matrices, i.e., for `alpha`>1 do not pass in generating matrices which are already interlaced. The Note for this class contains more info.	`1`
`msb`	`bool`	Flag for Most Significant Bit (MSB) vs Least Significant Bit (LSB) integer representations in generating matrices. If `msb=False` (LSB order), then integers in generating matrices will be bit-reversed.	`None`
`_verbose`	`bool`	If `True`, print linear matrix scrambling matrices.	`False`

Source code in qmcpy/discrete_distribution/digital_net_b2/digital_net_b2.py

def __init__(self,
             dimension = 1,
             replications = None,
             seed = None,
             randomize = 'LMS DS',
             generating_matrices = "joe_kuo.6.21201.txt",
             order = 'RADICAL INVERSE',
             t = 63,
             alpha = 1,
             msb = None,
             _verbose = False,
             # deprecated
             graycode = None,
             t_max = None,
             t_lms = None):
    r"""
    Args:
        dimension (Union[int,np.ndarray]): Dimension of the generator.

            - If an `int` is passed in, use generating vector components at indices 0,...,`dimension`-1.
            - If an `np.ndarray` is passed in, use generating vector components at these indices.

        replications (int): Number of independent randomizations of a pointset.
        seed (Union[None,int,np.random.SeedSeq): Seed the random number generator for reproducibility.
        randomize (str): Options are

            - `'LMS DS'`: Linear matrix scramble with digital shift.
            - `'LMS'`: Linear matrix scramble only.
            - `'DS'`: Digital shift only.
            - `'NUS'`: Nested uniform scrambling. Also known as Owen scrambling. 
            - `'FALSE'`: No randomization. In this case the first point will be the origin. 

        generating_matrices (Union[str,np.ndarray,int]: Specify the generating matrices.

            - A `str` should be the name (or path) of a file from the LDData repo at [https://github.com/QMCSoftware/LDData/tree/main/dnet](https://github.com/QMCSoftware/LDData/tree/main/dnet).
            - An `np.ndarray` of integers with shape $(d,m_\mathrm{max})$ or $(r,d,m_\mathrm{max})$ where $d$ is the number of dimensions, $r$ is the number of replications, and $2^{m_\mathrm{max}}$ is the maximum number of supported points. Setting `msb=False` will flip the bits of ints in the generating matrices.

        order (str): `'RADICAL INVERSE'`, or `'GRAY'` ordering. See the doctest example above.
        t (int): Number of bits in integer represetation of points *after* randomization. The number of bits in the generating matrices is inferred based on the largest value.
        alpha (int): Interlacing factor for higher order nets.  
            When `alpha`>1, interlacing is performed regardless of the generating matrices,  
            i.e., for `alpha`>1 do *not* pass in generating matrices which are already interlaced.  
            The Note for this class contains more info.  
        msb (bool): Flag for Most Significant Bit (MSB) vs Least Significant Bit (LSB) integer representations in generating matrices. If `msb=False` (LSB order), then integers in generating matrices will be bit-reversed. 
        _verbose (bool): If `True`, print linear matrix scrambling matrices. 
    """
    if graycode is not None:
        order = 'GRAY' if graycode else 'RADICAL INVERSE'
        warnings.warn("graycode argument deprecated, set order='GRAY' or order='RADICAL INVERSE' instead. Using order='%s'"%order,ParameterWarning)
    if t_lms is not None:
        t = t_lms
        warnings.warn("t_lms argument deprecated. Set t instead. Using t = %d"%t,ParameterWarning)
    if t_max is not None: 
        warnings.warn("t_max is deprecated as it can be inferred from the generating matrices. Set t to change the number of bits after randomization.",ParameterWarning)
    self.parameters = ['randomize','gen_mats_source','order','t','alpha','n_limit']
    self.input_generating_matrices = deepcopy(generating_matrices)
    self.input_t = deepcopy(t) 
    self.input_msb = deepcopy(msb)
    if isinstance(generating_matrices,str) and generating_matrices=="joe_kuo.6.21201.txt":
        self.gen_mats_source = generating_matrices
        if np.isscalar(dimension) and dimension<=1024 and alpha==1:
            gen_mats = np.load(dirname(abspath(__file__))+'/generating_matrices/joe_kuo.6.1024.npy')[None,:]
            d_limit = 1024
        else:
            gen_mats = np.load(dirname(abspath(__file__))+'/generating_matrices/joe_kuo.6.21201.npy')[None,:]
            d_limit = 21201
        msb = True
        n_limit = 4294967296
        self._t_curr = 32
        compat_shift = self._t_curr-t if self._t_curr>=t else 0
        if compat_shift>0: warnings.warn("Truncating ints in generating matrix to have t = %d bits."%t,ParameterWarning)
        gen_mats = gen_mats>>compat_shift
    elif isinstance(generating_matrices,str):
        self.gen_mats_source = generating_matrices
        assert generating_matrices[-4:]==".txt"
        local_root = dirname(abspath(__file__))+'/generating_matrices/'
        repos = DataSource()
        if repos.exists(local_root+generating_matrices):
            datafile = repos.open(local_root+generating_matrices)
        elif repos.exists("https://raw.githubusercontent.com/QMCSoftware/LDData/refs/heads/main/dnet/"+generating_matrices):
            datafile = repos.open("https://raw.githubusercontent.com/QMCSoftware/LDData/refs/heads/main/dnet/"+generating_matrices)
        elif repos.exists("https://raw.githubusercontent.com/QMCSoftware/LDData/refs/heads/main/"+generating_matrices):
            datafile = repos.open("https://raw.githubusercontent.com/QMCSoftware/LDData/refs/heads/main/"+generating_matrices)
        elif repos.exists("https://raw.githubusercontent.com/QMCSoftware/LDData/refs/heads/main/dnet/"+generating_matrices[7:]):
            datafile = repos.open("https://raw.githubusercontent.com/QMCSoftware/LDData/refs/heads/main/dnet/"+generating_matrices[7:])
        elif repos.exists("https://raw.githubusercontent.com/QMCSoftware/"+generating_matrices):
            datafile = repos.open("https://raw.githubusercontent.com/QMCSoftware/"+generating_matrices)
        elif repos.exists(generating_matrices):
            datafile = repos.open(generating_matrices)
        else:
            raise ParameterError("LDData path %s not found"%generating_matrices)
        contents = [line.rstrip('\n').strip() for line in datafile.readlines()]
        contents = [line.split("#",1)[0] for line in contents if line[0]!="#"]
        datafile.close()
        msb = True
        assert int(contents[0])==2, "DigitalNetB2 requires base=2 " # base 2
        d_limit = int(contents[1])
        n_limit = int(contents[2])
        self._t_curr = int(contents[3])
        compat_shift = self._t_curr-t if self._t_curr>=t else 0
        if compat_shift>0: warnings.warn("Truncating ints in generating matrix to have t = %d bits."%t,ParameterWarning)
        gen_mats = np.array([[int(v)>>compat_shift for v in line.split(' ')] for line in contents[4:]],dtype=np.uint64)[None,:]
    elif isinstance(generating_matrices,np.ndarray):
        self.gen_mats_source = "custom"
        assert generating_matrices.ndim==2 or generating_matrices.ndim==3
        gen_mats = generating_matrices[None,:,:] if generating_matrices.ndim==2 else generating_matrices
        assert isinstance(msb,bool), "when generating_matrices is a np.ndarray you must set either msb=True (for most significant bit ordering) or msb=False (for least significant bit ordering which will require a bit reversal)"
        gen_mat_max = gen_mats.max() 
        assert gen_mat_max>0, "generating matrix must have positive ints"
        self._t_curr = int(np.ceil(np.log2(gen_mat_max+1)))
        d_limit = gen_mats.shape[1]
        n_limit = int(2**(gen_mats.shape[2]))
    else:
        raise ParameterError("invalid generating_matrices, must be a string or np.ndarray.")
    super(DigitalNetB2,self).__init__(dimension,replications,seed,d_limit,n_limit)
    assert gen_mats.ndim==3 and gen_mats.shape[1]>=self.d and (gen_mats.shape[0]==1 or gen_mats.shape[0]==self.replications) and gen_mats.shape[2]>0, "invalid gen_mats.shape = %s"%str(gen_mats.shape)
    self.m_max = int(gen_mats.shape[-1])
    if isinstance(generating_matrices,np.ndarray) and msb:
        qmctoolscl.dnb2_gmat_lsb_to_msb(np.uint64(gen_mats.shape[0]),np.uint64(self.d),np.uint64(self.m_max),np.tile(np.uint64(self._t_curr),int(gen_mats.shape[0])),gen_mats,gen_mats,backend="c")
    self.order = str(order).upper().strip().replace("_"," ")
    if self.order=="GRAY CODE": self.order = "GRAY"
    if self.order=="NATURAL": self.order = "RADICAL INVERSE"
    assert self.order in ['RADICAL INVERSE','GRAY']
    assert isinstance(t,int) and t>0
    assert self._t_curr<=t<=64, "t must no more than 64 and no less than %d (the number of bits used to represent the generating matrices)"%(self._t_curr)
    assert isinstance(alpha,int) and alpha>0
    self.alpha = alpha
    if self.alpha>1:
        assert (self.dvec==np.arange(self.d)).all(), "digital interlacing requires dimension is an int"
        if self.m_max!=self._t_curr:
            warnings.warn("Digital interlacing is often performed on matrices with the number of columns (m_max = %d) equal to the number of bits in each int (%d), but this is not the case. Ensure you are NOT setting alpha>1 when generating matrices are already interlaced."%(self.m_max,self._t_curr),ParameterWarning)
    self._verbose = _verbose
    self.randomize = str(randomize).upper().strip().replace("_"," ")
    if self.randomize=="TRUE": self.randomize = "LMS DS"
    if self.randomize=="OWEN": self.randomize = "NUS"
    if self.randomize=="NONE": self.randomize = "FALSE"
    if self.randomize=="NO": self.randomize = "FALSE"
    assert self.randomize in ["LMS DS","LMS","DS","NUS","FALSE"]
    self.dtalpha = self.alpha*self.d
    if self.randomize=="FALSE":
        if self.alpha==1:
            self.gen_mats = gen_mats[:,self.dvec,:]
            self.t = self._t_curr
        else: 
            t_alpha = min(self.alpha*self._t_curr,t)
            gen_mat_ho = np.empty((gen_mats.shape[0],self.d,self.m_max),dtype=np.uint64)
            qmctoolscl.dnb2_interlace(np.uint64(gen_mats.shape[0]),np.uint64(self.d),np.uint64(self.m_max),np.uint64(self.dtalpha),np.uint64(self._t_curr),np.uint64(t_alpha),np.uint64(self.alpha),gen_mats[:,:self.dtalpha,:].copy(),gen_mat_ho,backend="c")
            self.gen_mats = gen_mat_ho
            self._t_curr = t_alpha
            self.t = self._t_curr
    elif self.randomize=="DS":
        if self.alpha==1:
            self.gen_mats = gen_mats[:,self.dvec,:]
            self.t = t
        else: 
            t_alpha = min(self.alpha*self._t_curr,t)
            gen_mat_ho = np.empty((gen_mats.shape[0],self.d,self.m_max),dtype=np.uint64)
            qmctoolscl.dnb2_interlace(np.uint64(gen_mats.shape[0]),np.uint64(self.d),np.uint64(self.m_max),np.uint64(self.dtalpha),np.uint64(self._t_curr),np.uint64(t_alpha),np.uint64(self.alpha),gen_mats[:,:self.dtalpha,:].copy(),gen_mat_ho,backend="c")
            self.gen_mats = gen_mat_ho
            self._t_curr = t_alpha
            self.t = t
        self.rshift = qmctoolscl.random_tbit_uint64s(self.rng,self.t,(self.replications,self.d))
    elif self.randomize in ["LMS","LMS DS"]:
        if self.alpha==1:
            gen_mat_lms = np.empty((self.replications,self.d,self.m_max),dtype=np.uint64)
            S = qmctoolscl.dnb2_get_linear_scramble_matrix(self.rng,np.uint64(self.replications),np.uint64(self.d),np.uint64(self._t_curr),np.uint64(t),np.uint64(self._verbose))
            qmctoolscl.dnb2_linear_matrix_scramble(np.uint64(self.replications),np.uint64(self.d),np.uint64(self.m_max),np.uint64(gen_mats.shape[0]),np.uint64(t),S,gen_mats[:,self.dvec,:].copy(),gen_mat_lms,backend="c")
            self.gen_mats = gen_mat_lms
            self._t_curr = t
            self.t = self._t_curr
        else:
            gen_mat_lms = np.empty((self.replications,self.dtalpha,self.m_max),dtype=np.uint64)
            S = qmctoolscl.dnb2_get_linear_scramble_matrix(self.rng,np.uint64(self.replications),np.uint64(self.dtalpha),np.uint64(self._t_curr),np.uint64(t),np.uint64(self._verbose))
            qmctoolscl.dnb2_linear_matrix_scramble(np.uint64(self.replications),np.uint64(self.dtalpha),np.uint64(self.m_max),np.uint64(gen_mats.shape[0]),np.uint64(t),S,gen_mats[:,:self.dtalpha,:].copy(),gen_mat_lms,backend="c")
            gen_mat_lms_ho = np.empty((self.replications,self.d,self.m_max),dtype=np.uint64)
            qmctoolscl.dnb2_interlace(np.uint64(self.replications),np.uint64(self.d),np.uint64(self.m_max),np.uint64(self.dtalpha),np.uint64(t),np.uint64(t),np.uint64(self.alpha),gen_mat_lms,gen_mat_lms_ho,backend="c")
            self.gen_mats = gen_mat_lms_ho
            self._t_curr = t
            self.t = self._t_curr
        if self.randomize=="LMS DS":
            self.rshift = qmctoolscl.random_tbit_uint64s(self.rng,self.t,(self.replications,self.d))
    elif self.randomize=="NUS":
        if alpha==1:
            new_seeds = self._base_seed.spawn(self.replications*self.d)
            self.rngs = np.array([np.random.Generator(np.random.SFC64(new_seeds[j])) for j in range(self.replications*self.d)]).reshape(self.replications,self.d)
            self.root_nodes = np.array([qmctoolscl.NUSNode_dnb2() for i in range(self.replications*self.d)]).reshape(self.replications,self.d)
            self.gen_mats = gen_mats[:,self.dvec,:].copy()
            self.t = t
        else:
            self.dtalpha = self.alpha*self.d
            new_seeds = self._base_seed.spawn(self.replications*self.dtalpha)
            self.rngs = np.array([np.random.Generator(np.random.SFC64(new_seeds[j])) for j in range(self.replications*self.dtalpha)]).reshape(self.replications,self.dtalpha)
            self.root_nodes = np.array([qmctoolscl.NUSNode_dnb2() for i in range(self.replications*self.dtalpha)]).reshape(self.replications,self.dtalpha)
            self.gen_mats = gen_mats[:,:self.dtalpha,:].copy()
            self.t = t
    else:
        raise ParameterError("self.randomize parsing error")
    self.gen_mats = np.ascontiguousarray(self.gen_mats)
    gen_mat_max = self.gen_mats.max() 
    assert gen_mat_max>0, "generating matrix must have positive ints"
    assert self._t_curr==int(np.ceil(np.log2(gen_mat_max+1)))
    assert 0<self._t_curr<=self.t<=64, "invalid 0 <= self._t_curr (%d) <= self.t (%d) <= 64"%(self._t_curr,self.t)
    if self.randomize=="FALSE": assert self.gen_mats.shape[0]==self.replications, "randomize='FALSE' but replications = %d does not equal the number of sets of generating matrices %d"%(self.replications,self.gen_mats.shape[0])

`Lattice`

Bases: AbstractLDDiscreteDistribution

Low discrepancy lattice sequence.

Note

Lattice sample sizes should be powers of \(2\) e.g. \(1\), \(2\), \(4\), \(8\), \(16\), \(\dots\).
The first point of an unrandomized lattice is the origin.

Examples:

>>> discrete_distrib = Lattice(2,seed=7)
>>> discrete_distrib(4)
array([[0.04386058, 0.58727432],
       [0.54386058, 0.08727432],
       [0.29386058, 0.33727432],
       [0.79386058, 0.83727432]])
>>> discrete_distrib(1) # first point in the sequence
array([[0.04386058, 0.58727432]])
>>> discrete_distrib
Lattice (AbstractLDDiscreteDistribution)
    d               2^(1)
    replications    1
    randomize       SHIFT
    gen_vec_source  kuo.lattice-33002-1024-1048576.9125.txt
    order           RADICAL INVERSE
    n_limit         2^(20)
    entropy         7

Replications of independent randomizations

>>> x = Lattice(3,seed=7,replications=2)(4)
>>> x.shape
(2, 4, 3)
>>> x
array([[[0.04386058, 0.58727432, 0.3691824 ],
        [0.54386058, 0.08727432, 0.8691824 ],
        [0.29386058, 0.33727432, 0.1191824 ],
        [0.79386058, 0.83727432, 0.6191824 ]],

       [[0.65212985, 0.69669968, 0.10605352],
        [0.15212985, 0.19669968, 0.60605352],
        [0.90212985, 0.44669968, 0.85605352],
        [0.40212985, 0.94669968, 0.35605352]]])

Different orderings (avoid warnings that the first point is the origin).

>>> Lattice(dimension=2,randomize=False,order='RADICAL INVERSE')(4,warn=False) 
array([[0.  , 0.  ],
       [0.5 , 0.5 ],
       [0.25, 0.75],
       [0.75, 0.25]])
>>> Lattice(dimension=2,randomize=False,order='GRAY')(4,warn=False)
array([[0.  , 0.  ],
       [0.5 , 0.5 ],
       [0.75, 0.25],
       [0.25, 0.75]])
>>> Lattice(dimension=2,randomize=False,order='LINEAR')(4,warn=False)
array([[0.  , 0.  ],
       [0.25, 0.75],
       [0.5 , 0.5 ],
       [0.75, 0.25]])

Generating vector from https://github.com/QMCSoftware/LDData/tree/main/lattice

>>> Lattice(dimension=3,randomize=False,generating_vector="mps.exod2_base2_m20_CKN.txt")(8,warn=False)
array([[0.   , 0.   , 0.   ],
       [0.5  , 0.5  , 0.5  ],
       [0.25 , 0.75 , 0.75 ],
       [0.75 , 0.25 , 0.25 ],
       [0.125, 0.375, 0.375],
       [0.625, 0.875, 0.875],
       [0.375, 0.125, 0.125],
       [0.875, 0.625, 0.625]])

Random generating vector supporting \(2^{25}\) points

>>> discrete_distrib = Lattice(3,generating_vector=25,seed=55,randomize=False)
>>> discrete_distrib.gen_vec
array([[       1, 11961679, 12107519]], dtype=uint64)
>>> discrete_distrib(4,warn=False)
array([[0.  , 0.  , 0.  ],
       [0.5 , 0.5 , 0.5 ],
       [0.25, 0.75, 0.75],
       [0.75, 0.25, 0.25]])

Two random generating vectors both supporting \(2^{25}\) points along with independent random shifts

>>> discrete_distrib = Lattice(3,seed=7,generating_vector=25,replications=2)
>>> discrete_distrib.gen_vec
array([[       1, 32809149,  1471719],
       [       1,   275319, 19705657]], dtype=uint64)
>>> discrete_distrib(4)
array([[[0.3691824 , 0.65212985, 0.69669968],
        [0.8691824 , 0.15212985, 0.19669968],
        [0.6191824 , 0.90212985, 0.44669968],
        [0.1191824 , 0.40212985, 0.94669968]],

       [[0.10605352, 0.63025643, 0.13630282],
        [0.60605352, 0.13025643, 0.63630282],
        [0.35605352, 0.38025643, 0.38630282],
        [0.85605352, 0.88025643, 0.88630282]]])

References

Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou.
GAIL: Guaranteed Automatic Integration Library (Version 2.3), MATLAB Software, 2019.
http://gailgithub.github.io/GAIL_Dev/.
F.Y. Kuo, D. Nuyens.
Application of quasi-Monte Carlo methods to elliptic PDEs with random diffusion coefficients - a survey of analysis and implementation.
Foundations of Computational Mathematics, 16(6):1631-1696, 2016.
https://link.springer.com/article/10.1007/s10208-016-9329-5.
D. Nuyens.
The Magic Point Shop of QMC point generators and generating vectors.
MATLAB and Python software, 2018.
https://people.cs.kuleuven.be/~dirk.nuyens/.
R. Cools, F.Y. Kuo, D. Nuyens.
Constructing embedded lattice rules for multivariate integration.
SIAM J. Sci. Comput., 28(6), 2162-2188.
P. L'Ecuyer, D. Munger.
LatticeBuilder: A General Software Tool for Constructing Rank-1 Lattice Rules.
ACM Transactions on Mathematical Software. 42. (2015).
10.1145/2754929.

Parameters:

Name	Type	Description	Default
`dimension`	`Union[int, ndarray]`	Dimension of the generator. If an `int` is passed in, use generating vector components at indices 0,...,`dimension`-1. If an `np.ndarray` is passed in, use generating vector components at these indices.	`1`
`replications`	`int`	Number of independent randomizations.	`None`
`seed`	`Union[None,int,np.random.SeedSeq`	Seed the random number generator for reproducibility.	`None`
`randomize`	`str`	Options are `'SHIFT'`: Random shift. `'FALSE'`: No randomization. In this case the first point will be the origin.	`'SHIFT'`
`generating_vector`	`Union[str, ndarray, int]`	Specify the generating vector. A `str` should be the name (or path) of a file from the LDData repo at https://github.com/QMCSoftware/LDData/tree/main/lattice. A `np.ndarray` of integers with shape \((d,)\) or \((r,d)\) where \(d\) is the number of dimensions and \(r\) is the number of replications. Must supply `m_max` where \(2^{m_\mathrm{max}}\) is the max number of supported samples. An `int`, call it \(M\), gives the random generating vector \((1,v_1,\dots,v_{d-1})^T\) where \(d\) is the dimension and \(v_i\) are randomly selected from \(\{3,5,\dots,2^M-1\}\) uniformly and independently. We require require \(1 < M < 27\).	`'kuo.lattice-33002-1024-1048576.9125.txt'`
`order`	`str`	`'LINEAR'`, `'RADICAL INVERSE'`, or `'GRAY'` ordering. See the doctest example above.	`'RADICAL INVERSE'`
`m_max`	`int`	\(2^{m_\mathrm{max}}\) is the maximum number of supported samples.	`None`

Source code in qmcpy/discrete_distribution/lattice/lattice.py

def __init__(self,
             dimension = 1,
             replications = None,
             seed = None,
             randomize = 'SHIFT',
             generating_vector = "kuo.lattice-33002-1024-1048576.9125.txt",
             order = 'RADICAL INVERSE',
             m_max = None):
    r"""
    Args:
        dimension (Union[int,np.ndarray]): Dimension of the generator.

            - If an `int` is passed in, use generating vector components at indices 0,...,`dimension`-1.
            - If an `np.ndarray` is passed in, use generating vector components at these indices.

        replications (int): Number of independent randomizations.
        seed (Union[None,int,np.random.SeedSeq): Seed the random number generator for reproducibility.
        randomize (str): Options are 

            - `'SHIFT'`: Random shift.
            - `'FALSE'`: No randomization. In this case the first point will be the origin. 

        generating_vector (Union[str,np.ndarray,int]: Specify the generating vector.

            - A `str` should be the name (or path) of a file from the LDData repo at [https://github.com/QMCSoftware/LDData/tree/main/lattice](https://github.com/QMCSoftware/LDData/tree/main/lattice).
            - A `np.ndarray` of integers with shape $(d,)$ or $(r,d)$ where $d$ is the number of dimensions and $r$ is the number of replications.  
                Must supply `m_max` where $2^{m_\mathrm{max}}$ is the max number of supported samples. 
            - An `int`, call it $M$, 
            gives the random generating vector $(1,v_1,\dots,v_{d-1})^T$ 
            where $d$ is the dimension and $v_i$ are randomly selected from $\{3,5,\dots,2^M-1\}$ uniformly and independently.  
            We require require $1 < M < 27$. 

        order (str): `'LINEAR'`, `'RADICAL INVERSE'`, or `'GRAY'` ordering. See the doctest example above.
        m_max (int): $2^{m_\mathrm{max}}$ is the maximum number of supported samples.
    """
    self.parameters = ['randomize','gen_vec_source','order','n_limit']
    self.input_generating_vector = deepcopy(generating_vector)
    self.input_m_max = deepcopy(m_max)
    if isinstance(generating_vector,str) and generating_vector=="kuo.lattice-33002-1024-1048576.9125.txt":
        self.gen_vec_source = generating_vector
        gen_vec = np.load(dirname(abspath(__file__))+'/generating_vectors/kuo.lattice-33002-1024-1048576.9125.npy')[None,:]
        d_limit = 9125
        n_limit = 1048576
    elif isinstance(generating_vector,str):
        self.gen_vec_source = generating_vector
        assert generating_vector[-4:]==".txt"
        local_root = dirname(abspath(__file__))+'/generating_vectors/'
        repos = DataSource()
        if repos.exists(local_root+generating_vector):
            datafile = repos.open(local_root+generating_vector)
        elif repos.exists("https://raw.githubusercontent.com/QMCSoftware/LDData/refs/heads/main/lattice/"+generating_vector):
            datafile = repos.open("https://raw.githubusercontent.com/QMCSoftware/LDData/refs/heads/main/lattice/"+generating_vector)
        elif repos.exists("https://raw.githubusercontent.com/QMCSoftware/LDData/refs/heads/main/"+generating_vector):
            datafile = repos.open("https://raw.githubusercontent.com/QMCSoftware/LDData/refs/heads/main/"+generating_vector)
        elif repos.exists("https://raw.githubusercontent.com/QMCSoftware/LDData/refs/heads/main/lattice/"+generating_vector[7:]):
            datafile = repos.open("https://raw.githubusercontent.com/QMCSoftware/LDData/refs/heads/main/lattice/"+generating_vector[7:])
        elif repos.exists("https://raw.githubusercontent.com/QMCSoftware/"+generating_vector):
            datafile = repos.open("https://raw.githubusercontent.com/QMCSoftware/"+generating_vector)
        elif repos.exists(generating_vector):
            datafile = repos.open(generating_vector)
        else:
            raise ParameterError("LDData path %s not found"%generating_vector)
        contents = [int(line.rstrip('\n').strip().split("#",1)[0]) for line in datafile.readlines() if line[0]!="#"]
        datafile.close()
        d_limit = int(contents[0])
        n_limit = int(contents[1])
        gen_vec = np.array(contents[2:],dtype=np.uint64)[None,:]
    elif isinstance(generating_vector,np.ndarray):
        self.gen_vec_source = "custom"
        gen_vec = generating_vector
        if m_max is None:
            raise ParameterError("m_max must be supplied when generating_vector is a np.ndarray")
        n_limit = int(2**m_max)
        d_limit = int(gen_vec.shape[-1])
    elif isinstance(generating_vector,int):
        assert 1<generating_vector<27, "int generating vector out of range"
        n_limit = 2**generating_vector
        assert isinstance(dimension,int), "random generating vector requires int dimension"
        d_limit = dimension
    else:
        raise ParameterError("invalid generating_vector, must be a string, numpy.ndarray, or int")
    super(Lattice,self).__init__(dimension,replications,seed,d_limit,n_limit)
    if isinstance(generating_vector,int):
        self.gen_vec_source = "random"
        m_max = int(np.log2(self.n_limit))
        gen_vec = np.hstack([np.ones((self.replications,1),dtype=np.uint64),2*self.rng.integers(1,2**(m_max-1),size=(self.replications,dimension-1),dtype=np.uint64)+1]).copy()
    assert isinstance(gen_vec,np.ndarray)
    gen_vec = np.atleast_2d(gen_vec) 
    assert gen_vec.ndim==2 and gen_vec.shape[1]>=self.d and (gen_vec.shape[0]==1 or gen_vec.shape[0]==self.replications), "invalid gen_vec.shape = %s"%str(gen_vec.shape)
    self.gen_vec = gen_vec[:,self.dvec].copy()
    self.order = str(order).upper().strip().replace("_"," ")
    if self.order=="GRAY CODE": self.order = "GRAY"
    if self.order=="NATURAL": self.order = "RADICAL INVERSE"
    assert self.order in ['LINEAR','RADICAL INVERSE','GRAY']
    self.randomize = str(randomize).upper()
    if self.randomize=="TRUE": self.randomize = "SHIFT"
    if self.randomize=="NONE": self.randomize = "FALSE"
    if self.randomize=="NO": self.randomize = "FALSE"
    assert self.randomize in ["SHIFT","FALSE"]
    if self.randomize=="SHIFT":
        self.shift = self.rng.uniform(size=(self.replications,self.d))
    if self.randomize=="FALSE": assert self.gen_vec.shape[0]==self.replications, "randomize='FALSE' but replications = %d does not equal the number of sets of generating vectors %d"%(self.replications,self.gen_vec.shape[0])

`Halton`

Bases: AbstractLDDiscreteDistribution

Low discrepancy Halton points.

Note

The first point of an unrandomized Halton sequence is the origin.
QRNG does not support multiple replications (independent randomizations).

Examples:

>>> discrete_distrib = Halton(2,seed=7)
>>> discrete_distrib(4)
array([[0.83790457, 0.89981478],
       [0.00986102, 0.4610941 ],
       [0.62236343, 0.02796307],
       [0.29427505, 0.79909098]])
>>> discrete_distrib
Halton (AbstractLDDiscreteDistribution)
    d               2^(1)
    replications    1
    randomize       LMS DP
    t               63
    n_limit         2^(32)
    entropy         7

Replications of independent randomizations

>>> x = Halton(3,seed=7,replications=2)(4)
>>> x.shape
(2, 4, 3)
>>> x
array([[[0.70988236, 0.18180876, 0.54073621],
        [0.38178158, 0.61168824, 0.64684354],
        [0.98597752, 0.70650871, 0.31479029],
        [0.15795399, 0.28162992, 0.98945647]],

       [[0.620398  , 0.57025403, 0.46336542],
        [0.44021889, 0.69926312, 0.60133428],
        [0.89132308, 0.12030255, 0.35715804],
        [0.04025218, 0.44304244, 0.10724799]]])

Unrandomized Halton

>>> Halton(2,randomize="FALSE",seed=7)(4,warn=False)
array([[0.        , 0.        ],
       [0.5       , 0.33333333],
       [0.25      , 0.66666667],
       [0.75      , 0.11111111]])

All randomizations

>>> Halton(2,randomize="LMS DP",seed=7)(4)
array([[0.83790457, 0.89981478],
       [0.00986102, 0.4610941 ],
       [0.62236343, 0.02796307],
       [0.29427505, 0.79909098]])
>>> Halton(2,randomize="LMS DS",seed=7)(4)
array([[0.82718745, 0.90603116],
       [0.0303368 , 0.44704107],
       [0.60182684, 0.03580544],
       [0.30505343, 0.78367016]])
>>> Halton(2,randomize="LMS",seed=7)(4,warn=False)
array([[0.        , 0.        ],
       [0.82822666, 0.92392942],
       [0.28838899, 0.46493682],
       [0.6165384 , 0.2493814 ]])
>>> Halton(2,randomize="DP",seed=7)(4)
array([[0.11593484, 0.99232505],
       [0.61593484, 0.65899172],
       [0.36593484, 0.32565839],
       [0.86593484, 0.77010283]])
>>> Halton(2,randomize="DS",seed=7)(4)
array([[0.56793849, 0.04063513],
       [0.06793849, 0.37396846],
       [0.81793849, 0.7073018 ],
       [0.31793849, 0.15174624]])
>>> Halton(2,randomize="NUS",seed=7)(4)
array([[0.141964  , 0.99285569],
       [0.65536579, 0.51938353],
       [0.46955206, 0.11342811],
       [0.78505432, 0.87032345]])
>>> Halton(2,randomize="QRNG",seed=7)(4)
array([[0.35362988, 0.38733489],
       [0.85362988, 0.72066823],
       [0.10362988, 0.05400156],
       [0.60362988, 0.498446  ]])

Replications of randomizations

>>> Halton(3,randomize="LMS DP",seed=7,replications=2)(4)
array([[[0.70988236, 0.18180876, 0.54073621],
        [0.38178158, 0.61168824, 0.64684354],
        [0.98597752, 0.70650871, 0.31479029],
        [0.15795399, 0.28162992, 0.98945647]],

       [[0.620398  , 0.57025403, 0.46336542],
        [0.44021889, 0.69926312, 0.60133428],
        [0.89132308, 0.12030255, 0.35715804],
        [0.04025218, 0.44304244, 0.10724799]]])
>>> Halton(3,randomize="LMS DS",seed=7,replications=2)(4)
array([[[4.57465163e-01, 5.75419751e-04, 7.47353067e-01],
        [6.29314800e-01, 9.24349881e-01, 8.47915779e-01],
        [2.37544271e-01, 4.63986168e-01, 1.78817056e-01],
        [9.09318567e-01, 2.48566227e-01, 3.17475640e-01]],

       [[6.04003127e-01, 9.92849835e-01, 4.21625151e-01],
        [4.57027115e-01, 1.97310094e-01, 2.43670150e-01],
        [8.76467351e-01, 4.22339232e-01, 1.05777101e-01],
        [5.46933622e-02, 7.79075280e-01, 9.29409300e-01]]])
>>> Halton(3,randomize="LMS",seed=7,replications=2)(4,warn=False)
array([[[0.        , 0.        , 0.        ],
        [0.82822666, 0.92392942, 0.34057871],
        [0.28838899, 0.46493682, 0.47954399],
        [0.6165384 , 0.2493814 , 0.77045601]],

       [[0.        , 0.        , 0.        ],
        [0.93115665, 0.57483093, 0.87170952],
        [0.48046642, 0.8122114 , 0.69381851],
        [0.58055977, 0.28006957, 0.55586147]]])
>>> Halton(3,randomize="DS",seed=7,replications=2)(4)
array([[[0.56793849, 0.04063513, 0.74276256],
        [0.06793849, 0.37396846, 0.94276256],
        [0.81793849, 0.7073018 , 0.14276256],
        [0.31793849, 0.15174624, 0.34276256]],

       [[0.98309816, 0.80260469, 0.17299622],
        [0.48309816, 0.13593802, 0.37299622],
        [0.73309816, 0.46927136, 0.57299622],
        [0.23309816, 0.9137158 , 0.77299622]]])
>>> Halton(3,randomize="DP",seed=7,replications=2)(4)
array([[[0.11593484, 0.99232505, 0.6010751 ],
        [0.61593484, 0.65899172, 0.0010751 ],
        [0.36593484, 0.32565839, 0.4010751 ],
        [0.86593484, 0.77010283, 0.8010751 ]],

       [[0.26543198, 0.12273092, 0.20202896],
        [0.76543198, 0.45606426, 0.60202896],
        [0.01543198, 0.78939759, 0.40202896],
        [0.51543198, 0.23384203, 0.00202896]]])
>>> Halton(3,randomize="NUS",seed=7,replications=2)(4)
array([[[0.141964  , 0.99285569, 0.77722918],
        [0.65536579, 0.51938353, 0.22797442],
        [0.46955206, 0.11342811, 0.9975298 ],
        [0.78505432, 0.87032345, 0.57696123]],

       [[0.04813634, 0.16158904, 0.56038465],
        [0.89364888, 0.33578478, 0.36145822],
        [0.34111023, 0.84596814, 0.0292313 ],
        [0.71866903, 0.23852281, 0.80431142]]])

References:

Marius Hofert and Christiane Lemieux.
qrng: (Randomized) Quasi-Random Number Generators.
R package version 0.0-7. (2019).
https://CRAN.R-project.org/package=qrng.
A. B. Owen.
A randomized Halton algorithm in R.
arXiv:1706.02808 [stat.CO]. 2017.
A. B. Owen and Z. Pan.
Gain coefficients for scrambled Halton points.
arXiv:2308.08035 [stat.CO]. 2023.

Parameters:

Name	Type	Description	Default
`dimension`	`Union[int, ndarray]`	Dimension of the generator. If an `int` is passed in, use generating vector components at indices 0,...,`dimension`-1. If an `np.ndarray` is passed in, use generating vector components at these indices.	`1`
`replications`	`int`	Number of independent randomizations of a pointset.	`None`
`seed`	`Union[None,int,np.random.SeedSeq`	Seed the random number generator for reproducibility.	`None`
`randomize`	`str`	Options are `'LMS DP'`: Linear matrix scramble with digital permutation. `'LMS DS'`: Linear matrix scramble with digital shift. `'LMS'`: Linear matrix scramble only. `'DP'`: Digital permutation scramble only. `'DS'`: Digital shift only. `'NUS'`: Nested uniform scrambling. `'QRNG'`: Deterministic permutation scramble and random digital shift from QRNG [1] (with `generalize=True`). Does not support replications>1. `None`: No randomization. In this case the first point will be the origin.	`'LMS DP'`
`t`	`int`	Number of bits in integer represetation of points after randomization. The number of bits in the generating matrices is inferred based on the largest value.	`63`
`n_lim`	`int`	Maximum number of compatible points, determines the number of rows in the generating matrices.	`2 ** 32`

Source code in qmcpy/discrete_distribution/halton.py

def __init__(self,
             dimension = 1,
             replications = None,
             seed = None, 
             randomize = 'LMS DP',
             t = 63,
             n_lim = 2**32,
             # deprecated
             t_lms = None):
    r"""
    Args:
        dimension (Union[int,np.ndarray]): Dimension of the generator.

            - If an `int` is passed in, use generating vector components at indices 0,...,`dimension`-1.
            - If an `np.ndarray` is passed in, use generating vector components at these indices.

        replications (int): Number of independent randomizations of a pointset.
        seed (Union[None,int,np.random.SeedSeq): Seed the random number generator for reproducibility.
        randomize (str): Options are

            - `'LMS DP'`: Linear matrix scramble with digital permutation.
            - `'LMS DS'`: Linear matrix scramble with digital shift.
            - `'LMS'`: Linear matrix scramble only.
            - `'DP'`: Digital permutation scramble only.
            - `'DS'`: Digital shift only.
            - `'NUS'`: Nested uniform scrambling.
            - `'QRNG'`: Deterministic permutation scramble and random digital shift from QRNG [1] (with `generalize=True`). Does *not* support replications>1.
            - `None`: No randomization. In this case the first point will be the origin. 
        t (int): Number of bits in integer represetation of points *after* randomization. The number of bits in the generating matrices is inferred based on the largest value.
        n_lim (int): Maximum number of compatible points, determines the number of rows in the generating matrices. 
    """
    if t_lms is not None:
        t = t_lms
        warnings.warn("t_lms argument deprecated. Set t instead. Using t = %d"%t,ParameterWarning)
    self.parameters = ['randomize','t','n_limit']
    self.all_primes = np.array([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919],dtype=np.uint64)
    d_limit = len(self.all_primes)
    self.input_t = deepcopy(t) 
    super(Halton,self).__init__(dimension,replications,seed,d_limit,n_lim)
    self.randomize = str(randomize).upper().strip().replace("_"," ")
    if self.randomize=="TRUE": self.randomize = "LMS DP"
    if self.randomize=="LMS PERM": self.randomize = "LMS DP"
    if self.randomize=="PERM": self.randomize = "DP"
    if self.randomize=="OWEN": self.randomize = "NUS"
    if self.randomize=="NONE": self.randomize = "FALSE"
    if self.randomize=="NO": self.randomize = "FALSE"
    assert self.randomize in ["LMS DP","LMS DS","LMS","DP","DS","NUS","QRNG","FALSE"]
    if self.randomize=="QRNG":
        from ._c_lib import _load_c_lib
        assert self.replications==1, "QRNG requires replications=1"
        self.randu_d_32 = self.rng.uniform(size=(self.d,32))
        _c_lib = _load_c_lib()
        import ctypes
        self.halton_cf_qrng = _c_lib.halton_qrng
        self.halton_cf_qrng.argtypes = [
            ctypes.c_int,  # n
            ctypes.c_int,  # d
            ctypes.c_int, # n0
            ctypes.c_int, # generalized
            np.ctypeslib.ndpointer(ctypes.c_double, flags='C_CONTIGUOUS'), # res
            np.ctypeslib.ndpointer(ctypes.c_double, flags='C_CONTIGUOUS'), # randu_d_32
            np.ctypeslib.ndpointer(ctypes.c_int, flags='C_CONTIGUOUS')]  # dvec
        self.halton_cf_qrng.restype = None
    self.primes = self.all_primes[self.dvec]
    self.m_max = int(np.ceil(np.log(self.n_limit)/np.log(self.primes.min())))
    self._t_curr = self.m_max
    self.t = self.m_max if self.m_max>t else t
    self.C = qmctoolscl.gdn_get_halton_generating_matrix(np.uint64(1),np.uint64(self.d),np.uint64(self._t_curr))
    if "LMS" in self.randomize:
        S = qmctoolscl.gdn_get_linear_scramble_matrix(self.rng,np.uint64(self.replications),np.uint64(self.d),np.uint64(self._t_curr),np.uint64(self.t),np.uint64(1),self.primes)
        C_lms = np.empty((self.replications,self.d,self.m_max,self.t),dtype=np.uint64)
        qmctoolscl.gdn_linear_matrix_scramble(np.uint64(self.replications),np.uint64(self.d),np.uint64(self.m_max),np.uint64(1),np.uint64(1),np.uint64(self._t_curr),np.uint64(self.t),self.primes,S,self.C,C_lms,backend="c")
        self.C = C_lms
        self._t_curr = self.t
    if "DP" in self.randomize:
        self.perms = qmctoolscl.gdn_get_digital_permutations(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(1),self.primes)
    if "DS" in self.randomize:
        self.rshift = qmctoolscl.gdn_get_digital_shifts(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(1),self.primes)
    if "NUS" in self.randomize:
        new_seeds = self._base_seed.spawn(self.replications*self.d)
        self.rngs = np.array([np.random.Generator(np.random.SFC64(new_seeds[j])) for j in range(self.replications*self.d)]).reshape(self.replications,self.d)
        self.root_nodes = np.array([qmctoolscl.NUSNode_gdn() for i in range(self.replications*self.d)]).reshape(self.replications,self.d)
    assert self.C.ndim==4 and (self.C.shape[0]==1 or self.C.shape[0]==self.replications) and self.C.shape[1]==self.d and self.C.shape[2]==self.m_max and self.C.shape[3]==self._t_curr
    assert 0<self._t_curr<=self.t<=64
    if self.randomize=="FALSE": assert self.C.shape[0]==self.replications, "randomize='FALSE' but replications = %d does not equal the number of sets of generating vectors %d"%(self.replications,self.C.shape[0])

`IIDStdUniform`

Bases: AbstractIIDDiscreteDistribution

IID standard uniform points, a wrapper around numpy.random.rand.

Note

Unlike low discrepancy sequence, calling an IIDStdUniform instance gives new samples every time,
e.g., running the first doctest below with dd = Lattice(dimension=2) would give the same 4 points in both calls,
but since we are using an IIDStdUniform instance it gives different points every call.

Examples:

>>> discrete_distrib = IIDStdUniform(dimension=2,seed=7)
>>> discrete_distrib(4)
array([[0.04386058, 0.58727432],
       [0.3691824 , 0.65212985],
       [0.69669968, 0.10605352],
       [0.63025643, 0.13630282]])
>>> discrete_distrib(4) # gives new samples every time
array([[0.5968363 , 0.0576251 ],
       [0.2028797 , 0.22909681],
       [0.1366783 , 0.75220658],
       [0.84501765, 0.56269008]])
>>> discrete_distrib
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               2^(1)
    replications    1
    entropy         7

Replications (implemented for API consistency)

>>> x = IIDStdUniform(dimension=3,replications=2,seed=7)(4)
>>> x.shape
(2, 4, 3)
>>> x
array([[[0.04386058, 0.58727432, 0.3691824 ],
        [0.65212985, 0.69669968, 0.10605352],
        [0.63025643, 0.13630282, 0.5968363 ],
        [0.0576251 , 0.2028797 , 0.22909681]],

       [[0.1366783 , 0.75220658, 0.84501765],
        [0.56269008, 0.04826852, 0.71308655],
        [0.80983568, 0.85383675, 0.80475135],
        [0.6171181 , 0.1239209 , 0.16809479]]])

Parameters:

Name	Type	Description	Default
`dimension`	`int`	Dimension of the samples.	`1`
`replications`	`Union[None, int]`	Number of randomizations. This is implemented only for API consistency. Equivalent to reshaping samples.	`None`
`seed`	`Union[None,int,np.random.SeedSeq`	Seed the random number generator for reproducibility.	`None`

Source code in qmcpy/discrete_distribution/iid_std_uniform.py

def __init__(self,
             dimension = 1, 
             replications = None, 
             seed = None):
    r"""
    Args:
        dimension (int): Dimension of the samples.
        replications (Union[None,int]): Number of randomizations. This is implemented only for API consistency. Equivalent to reshaping samples. 
        seed (Union[None,int,np.random.SeedSeq): Seed the random number generator for reproducibility.
    """
    super(IIDStdUniform,self).__init__(int(dimension),replications,seed,d_limit=np.inf,n_limit=np.inf)
    if not (self.dvec==np.arange(self.d)).all():
        warnings.warn("IIDStdUniform does not accomodate dvec",ParameterWarning)

Discrete Distributions

UML Overview

AbstractDiscreteDistribution

__call__

spawn

DigitalNetB2

Lattice

Halton

IIDStdUniform

UML Specific

`AbstractDiscreteDistribution`

call

`DigitalNetB2`

`Lattice`

`Halton`

`IIDStdUniform`