Skip to content

Discrete Distributions

UML Overview

AbstractDiscreteDistribution

Bases: object

Source code in qmcpy/discrete_distribution/abstract_discrete_distribution.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
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 
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
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 
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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:

  1. 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.

  2. 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.

  3. 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.

  4. 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/.

  5. R. Cools, F.Y. Kuo, D. Nuyens. Constructing embedded lattice rules for multivariate integration. SIAM J. Sci. Comput., 28(6), 2162-2188.

  6. 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).

  7. Sobol, Ilya & Asotsky, Danil & Kreinin, Alexander & Kucherenko, Sergei. (2011). Construction and Comparison of High-Dimensional Sobol' Generators. Wilmott. 2011. 10.1002/wilm.10056.

  8. 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 
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
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 (not 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

  1. 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/.

  2. 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.

  3. 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/.

  4. R. Cools, F.Y. Kuo, D. Nuyens. Constructing embedded lattice rules for multivariate integration. SIAM J. Sci. Comput., 28(6), 2162-2188.

  5. 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 
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
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: DigitalNetAnyBases

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:

  1. 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.

  2. A. B. Owen.
    A randomized Halton algorithm in R.
    arXiv:1706.02808 [stat.CO]. 2017.

  3. A. B. Owen and Z. Pan.
    Gain coefficients for scrambled Halton points.
    arXiv:2308.08035 [stat.CO]. 2023.

Source code in qmcpy/discrete_distribution/digital_net_any_bases/digital_net_any_bases.py 
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
def __init__(self,
             dimension = 1,
             replications = None,
             seed = None, 
             randomize = 'LMS DP',
             bases_generating_matrices = None,
             t = None,
             alpha = 1,
             n_lim = 2**32,
             warn = True):
    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. 
        bases_generating_matrices (Union[str, tuple]: Specify the bases and the generating matrices.

            - `"HALTON"` will use Halton generating matrices.
            - `"FAURE"` will use Faure generating matrices .
            - `bases,generating_matrices` requires 

                - `bases` is an `np.ndarray` of integers with shape $(,d)$ or $(r,d)$ where $d$ is the number of dimensions and $r$ is the number of replications.
                - `generating_matrices` is an `np.ndarray` of integers with shape $(d,m_\mathrm{max},t_\mathrm{max})$ or $(r,d,m_\mathrm{max},t_\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.

        t (int): Number of digits *after* randomization. The number of digits in the generating matrices is inferred.
        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.  
        n_lim (int): Maximum number of compatible points, determines the number of rows in the generating matrices. 
        warn (bool): If `False`, suppress warnings in construction 
    """
    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)
    if bases_generating_matrices is None:
        if self.DEFAULT_GENERATING_MATRICES=="HALTON":
            self.type_bases_generating_matrices = "HALTON"
            d_limit = len(self.all_primes)
        elif self.DEFAULT_GENERATING_MATRICES=="FAURE":
            self.type_bases_generating_matrices = "FAURE"
            d_limit = int(self.all_primes[-1])
        else:
            raise ParameterError("must supply bases_generating_matrices")
    else:
        self.type_bases_generating_matrices = "CUSTOM"
        assert len(bases_generating_matrices)==2
        bases,generating_matrices = bases_generating_matrices
        assert isinstance(generating_matrices,np.ndarray)
        assert generating_matrices.ndim==3 or generating_matrices.ndim==4 
        d_limit = generating_matrices.shape[1]
        if np.isscalar(bases):
            assert bases>0
            assert bases%1==0 
            bases = int(bases)*np.ones(d_limit,dtype=int)
        assert bases.ndim==1 or bases.ndim==2
    self.input_t = deepcopy(t) 
    self.input_bases_generating_matrices = deepcopy(bases_generating_matrices)
    super(DigitalNetAnyBases,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":
        assert self.type_bases_generating_matrices=="HALTON", "QRNG randomization is only applicable for the Halton generator."
        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.alpha = alpha
    assert self.alpha>=1
    assert self.alpha%1==0
    if self.alpha>1:
        assert (self.dvec==np.arange(self.d)).all(), "digital interlacing requires dimension is an int"
    self.dtalpha = self.alpha*self.d 
    if self.type_bases_generating_matrices=="HALTON":
        self.bases = self.all_primes[self.dvec][None,:]
        self.m_max = int(np.ceil(np.log(self.n_limit)/np.log(self.bases.min())))
        self._t_curr = self.m_max
        if t is None: t = int(np.ceil(-np.log(2**(-63))/np.log(self.bases.min())))
        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))
    elif self.type_bases_generating_matrices=="FAURE":
        assert (self.dvec==np.arange(self.d)).all(), "Faure requires dimension is an int"
        p = self.all_primes[np.argmax(self.all_primes>=self.d)]
        self.bases = p*np.ones((1,self.dtalpha),dtype=np.uint64)
        self.m_max = int(np.ceil(np.log(self.n_limit)/np.log(p)))
        if t is None: t = int(np.ceil(-np.log(2**(-63))/np.log(p)))
        self._t_curr = self.m_max
        self.t = self.m_max if self.m_max>t else t
        self.C = np.ones((self.dtalpha,1,1),dtype=np.uint64)*np.eye(self._t_curr,dtype=np.uint64)
        if self.dtalpha>1:
            for a in range(self._t_curr):
                for b in range(a+1):
                    self.C[1,a,b] = comb(a,b)%p
        if self.dtalpha>2:
            for k in range(2,self.dtalpha):
                for a in range(self._t_curr):
                    for b in range(a+1):
                        self.C[k,a,b] = (int(self.C[1,a,b])*((k**(a-b))%p))%p
        self.C = self.C[None,:,:,:]
        self.faure_base = p
    else:
        self.bases = bases.astype(np.uint64)
        if self.bases.ndim==1: self.bases = self.bases[None,:]
        assert self.bases.shape[1]>=self.dtalpha
        if self.alpha==1:
            self.bases = self.bases[:,self.dvec]
        else:
            self.bases = self.bases[:,:self.dtalpha]
        self.C = generating_matrices.astype(np.uint64)
        if self.C.ndim==3: self.C = self.C[None,:,:,:]
        assert self.C.shape[1]>=self.dtalpha
        if self.alpha==1:
            self.C = self.C[:,self.dvec,:,:]
        else:
            self.C = self.C[:,:self.dtalpha,:,:]
        self.m_max,self._t_curr = self.C.shape[-2:]
        if t is None: t = int(np.ceil(-np.log(2**(-63))/np.log(self.bases.min())))
        self.t = self.m_max if self.m_max>t else t
        assert (0<=self.C).all()
        assert (self.C<self.bases[:,:,None,None]).all()
        self.C = np.ascontiguousarray(self.C) 
        self.bases = np.ascontiguousarray(self.bases)
    if self.alpha>1:
        assert (self.bases==self.bases[0,0]).all(), "alpha>1 performs digital interlacing which requires the same base across dimensions and replications."
        if warn and self.m_max!=self._t_curr:
            warnings.warn("Digital interlacing is often performed on generating matrices with the number of columns (m_max = %d) equal to the number of rows (_t_curr = %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)
    assert self.bases.ndim==2
    assert self.bases.shape[-1]==self.dtalpha
    assert self.bases.shape[0]==1 or self.bases.shape[0]==self.replications
    assert self.C.ndim==4
    assert self.C.shape[-3:]==(self.dtalpha,self.m_max,self._t_curr)
    assert self.C.shape[0]==1 or self.C.shape[0]==self.replications
    r_b = self.bases.shape[0]
    r_C = self.C.shape[0]
    if self.randomize=="FALSE":
        if self.alpha>1:
            t_alpha = min(self.t,self._t_curr*self.alpha)
            C_ho = np.empty((self.replications,self.d,self.m_max,t_alpha),dtype=np.uint64)
            qmctoolscl.gdn_interlace(np.uint64(self.replications),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),self.C,C_ho)
            self.C = C_ho
            self._t_curr = t_alpha
            self.t = t_alpha
            self.bases = self.bases[:,:self.d]
    elif self.randomize=="DP":
        self.perms = qmctoolscl.gdn_get_digital_permutations(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(r_b),self.bases)
    elif self.randomize=="DS":
        self.rshift = qmctoolscl.gdn_get_digital_shifts(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(r_b),self.bases)
    elif self.randomize in ["LMS","LMS DS","LMS DP"]:
        if self.alpha==1:
            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(r_b),self.bases)
            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(r_C),np.uint64(r_b),np.uint64(self._t_curr),np.uint64(self.t),self.bases,S,self.C,C_lms,backend="c")
            self.C = C_lms
            self._t_curr = self.t
        else:
            t_dig = np.ceil(max(self.t/self.alpha,self._t_curr))
            S = qmctoolscl.gdn_get_linear_scramble_matrix(self.rng,np.uint64(self.replications),np.uint64(self.dtalpha),np.uint64(self._t_curr),np.uint64(t_dig),np.uint64(r_b),self.bases)
            C_lms = np.empty((self.replications,self.dtalpha,self.m_max,self.t),dtype=np.uint64)
            qmctoolscl.gdn_linear_matrix_scramble(np.uint64(self.replications),np.uint64(self.dtalpha),np.uint64(self.m_max),np.uint64(r_C),np.uint64(r_b),np.uint64(self._t_curr),np.uint64(t_dig),self.bases,S,self.C,C_lms,backend="c")
            C_lms_ho = np.empty((self.replications,self.d,self.m_max,self.t),dtype=np.uint64)
            qmctoolscl.gdn_interlace(np.uint64(self.replications),np.uint64(self.d),np.uint64(self.m_max),np.uint64(self.dtalpha),np.uint64(t_dig),np.uint64(self.t),np.uint64(self.alpha),C_lms,C_lms_ho)
            self.C = C_lms_ho
            self._t_curr = self.t
            self.t = self._t_curr
            self.bases = self.bases[:,:self.d]
        if self.randomize=="LMS DP":
            self.perms = qmctoolscl.gdn_get_digital_permutations(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(r_b),self.bases)
        elif self.randomize=="LMS DS":
            self.rshift = qmctoolscl.gdn_get_digital_shifts(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(r_b),self.bases)
    elif self.randomize=="NUS":
        if self.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_gdn() for i in range(self.replications*self.d)]).reshape(self.replications,self.d)
        else:
            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_gdn() for i in range(self.replications*self.dtalpha)]).reshape(self.replications,self.dtalpha)
    assert self.C.ndim==4 and (self.C.shape[0]==1 or self.C.shape[0]==self.replications) and self.C.shape[1]==(self.dtalpha if self.randomize=="NUS" else self.d) and self.C.shape[2]==self.m_max and self.C.shape[3]==self._t_curr
    assert self.bases.ndim==2 and (self.bases.shape[0]==1 or self.bases.shape[0]==self.replications) and self.bases.shape[1]==(self.dtalpha if self.randomize=="NUS" else self.d)
    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])
    if warn and (self.bases==2).all():
        warnings.warn("It is more efficient to use DigitalNetB2 instead of DigitalNetAnyBases when all bases are 2")
    self.warn = warn

Faure

Bases: DigitalNetAnyBases

Low discrepancy Faure points.

Note
  • The first point of an unrandomized Faure sequence is the origin.

Examples:

>>> discrete_distrib = Faure(4,seed=7)
>>> discrete_distrib(25)
array([[0.33534484, 0.26765076, 0.49135906, 0.6052543 ],
       [0.45194823, 0.53502655, 0.93195248, 0.5024511 ],
       [0.69367212, 0.09800441, 0.77275989, 0.35304609],
       [0.0086832 , 0.96517227, 0.01157255, 0.85055011],
       [0.97041531, 0.63254601, 0.25235643, 0.18869841],
       [0.20862319, 0.45844039, 0.64733082, 0.97109116],
       [0.57003274, 0.12560935, 0.08772865, 0.02771165],
       [0.73495613, 0.91292101, 0.32692686, 0.68483526],
       [0.05156003, 0.78796242, 0.56754229, 0.54342666],
       [0.89329109, 0.21506683, 0.80647015, 0.27293526],
       [0.25277941, 0.07296122, 0.23030907, 0.46332812],
       [0.4941915 , 0.94839898, 0.46928813, 0.39348452],
       [0.60952104, 0.69550381, 0.90983954, 0.81009581],
       [0.17092444, 0.37887902, 0.74911456, 0.10721583],
       [0.93584834, 0.40585698, 0.18944839, 0.72587572],
       [0.36989645, 0.85447136, 0.97641039, 0.14826857],
       [0.53450751, 0.7376398 , 0.61725827, 0.64572148],
       [0.65143961, 0.28501343, 0.05601974, 0.58393268],
       [0.09284913, 0.59218182, 0.29689116, 0.23247628],
       [0.80817253, 0.02749358, 0.53742056, 0.88960099],
       [0.29406032, 0.64634796, 0.15459067, 0.31205785],
       [0.40906423, 0.35338995, 0.39383649, 0.93066193],
       [0.77079531, 0.50849426, 0.43438996, 0.06810554],
       [0.13540741, 0.17593145, 0.87500744, 0.76631324],
       [0.85233694, 0.81903638, 0.71417647, 0.42286144]])
>>> discrete_distrib
Faure (AbstractLDDiscreteDistribution)
    d               2^(2)
    replications    1
    randomize       LMS DP
    t               28
    n_limit         2^(32)
    entropy         7

Replications of independent randomizations 

>>> x = Faure(3,seed=7,replications=2)(9)
>>> x.shape
(2, 9, 3)
>>> x
array([[[0.46995809, 0.81347921, 0.84921511],
        [0.22573669, 0.29702269, 0.03044283],
        [0.68125682, 0.38949622, 0.65611464],
        [0.58494536, 0.14367836, 0.43908813],
        [0.00739063, 0.63809438, 0.95363225],
        [0.79609189, 0.71822727, 0.1472069 ],
        [0.36690976, 0.47240941, 0.2470566 ],
        [0.12284101, 0.96819212, 0.55174327],
        [0.91154206, 0.05940035, 0.73295342]],

       [[0.44633162, 0.60684282, 0.96393795],
        [0.13904548, 0.73709456, 0.37711571],
        [0.79070871, 0.15560342, 0.12328662],
        [0.40389155, 0.97863405, 0.0114599 ],
        [0.05951757, 0.10889081, 0.85272281],
        [0.71514612, 0.52312704, 0.59877454],
        [0.64931396, 0.23505158, 0.48710589],
        [0.30097968, 0.36957094, 0.23358374],
        [0.99369356, 0.78380717, 0.74090153]]])

Unrandomized Faure 

>>> Faure(4,randomize="FALSE",seed=7)(25,warn=False)
array([[0.  , 0.  , 0.  , 0.  ],
       [0.2 , 0.2 , 0.2 , 0.2 ],
       [0.4 , 0.4 , 0.4 , 0.4 ],
       [0.6 , 0.6 , 0.6 , 0.6 ],
       [0.8 , 0.8 , 0.8 , 0.8 ],
       [0.04, 0.24, 0.44, 0.64],
       [0.24, 0.44, 0.64, 0.84],
       [0.44, 0.64, 0.84, 0.04],
       [0.64, 0.84, 0.04, 0.24],
       [0.84, 0.04, 0.24, 0.44],
       [0.08, 0.48, 0.88, 0.28],
       [0.28, 0.68, 0.08, 0.48],
       [0.48, 0.88, 0.28, 0.68],
       [0.68, 0.08, 0.48, 0.88],
       [0.88, 0.28, 0.68, 0.08],
       [0.12, 0.72, 0.32, 0.92],
       [0.32, 0.92, 0.52, 0.12],
       [0.52, 0.12, 0.72, 0.32],
       [0.72, 0.32, 0.92, 0.52],
       [0.92, 0.52, 0.12, 0.72],
       [0.16, 0.96, 0.76, 0.56],
       [0.36, 0.16, 0.96, 0.76],
       [0.56, 0.36, 0.16, 0.96],
       [0.76, 0.56, 0.36, 0.16],
       [0.96, 0.76, 0.56, 0.36]])

All randomizations 

>>> Faure(3,randomize="LMS DP",seed=7)(9)
array([[0.60869072, 0.76096155, 0.79807281],
       [0.03946861, 0.56443723, 0.39462162],
       [0.84046405, 0.17459934, 0.31263449],
       [0.39080182, 0.42253604, 0.09713822],
       [0.15491328, 0.22747922, 0.90356545],
       [0.95606114, 0.84998662, 0.50057213],
       [0.50599254, 0.09792123, 0.60600892],
       [0.26995157, 0.88915344, 0.20257439],
       [0.73776586, 0.51292533, 0.68846996]])
>>> Faure(3,randomize="LMS DS",seed=7)(9)
array([[0.09363194, 0.63822128, 0.52606134],
       [0.66021198, 0.02837507, 0.7155204 ],
       [0.85642224, 0.83340554, 0.22549083],
       [0.31989976, 0.30119144, 0.00893035],
       [0.55299403, 0.70247804, 0.63048602],
       [0.7492045 , 0.49633053, 0.81948782],
       [0.20034171, 0.96431985, 0.92797189],
       [0.43358859, 0.3653974 , 0.11741381],
       [0.9629798 , 0.17028087, 0.41752644]])
>>> Faure(3,randomize="LMS",seed=7)(9,warn=False)
array([[0.        , 0.        , 0.        ],
       [0.5706961 , 0.84756013, 0.63390413],
       [0.80394278, 0.65243798, 0.82336297],
       [0.22678173, 0.67547662, 0.93144564],
       [0.46399184, 0.50916041, 0.12088811],
       [0.69723852, 0.31536484, 0.40911072],
       [0.11956371, 0.33845241, 0.52953528],
       [0.35692624, 0.18443172, 0.70664719],
       [0.92335406, 0.97711587, 0.22942236]])
>>> Faure(3,randomize="DP",seed=7)(9)
array([[0.62331102, 0.65595286, 0.5208953 ],
       [0.28997769, 0.98928619, 0.18756197],
       [0.95664436, 0.32261952, 0.85422863],
       [0.51219991, 0.87817508, 0.96533975],
       [0.17886658, 0.21150841, 0.63200641],
       [0.84553325, 0.54484175, 0.29867308],
       [0.4010888 , 0.1003973 , 0.07645086],
       [0.06775547, 0.43373064, 0.74311752],
       [0.73442213, 0.76706397, 0.40978419]])
>>> Faure(3,randomize="DS",seed=7)(9)
array([[0.68058188, 0.90130449, 0.36582045],
       [0.01391521, 0.23463782, 0.69915378],
       [0.34724854, 0.56797116, 0.03248711],
       [0.79169299, 0.0124156 , 0.14359822],
       [0.12502632, 0.34574893, 0.47693156],
       [0.45835966, 0.67908227, 0.81026489],
       [0.9028041 , 0.45686005, 0.921376  ],
       [0.23613743, 0.79019338, 0.25470933],
       [0.56947077, 0.12352671, 0.58804267]])
>>> Faure(3,randomize="NUS",seed=7)(9)
array([[0.077534  , 0.99285569, 0.64186774],
       [0.82225806, 0.4963679 , 0.86453228],
       [0.48498632, 0.0833722 , 0.06011164],
       [0.15691802, 0.60236449, 0.11226224],
       [0.95445636, 0.2796012 , 0.54752025],
       [0.61890325, 0.84183901, 0.70654863],
       [0.25089638, 0.17805972, 0.95988146],
       [0.68344029, 0.77065782, 0.26676153],
       [0.4322891 , 0.40799837, 0.34911626]])

Replications of randomizations 

>>> Faure(3,randomize="LMS DP",seed=7,replications=2)(9)
array([[[0.46995809, 0.81347921, 0.84921511],
        [0.22573669, 0.29702269, 0.03044283],
        [0.68125682, 0.38949622, 0.65611464],
        [0.58494536, 0.14367836, 0.43908813],
        [0.00739063, 0.63809438, 0.95363225],
        [0.79609189, 0.71822727, 0.1472069 ],
        [0.36690976, 0.47240941, 0.2470566 ],
        [0.12284101, 0.96819212, 0.55174327],
        [0.91154206, 0.05940035, 0.73295342]],

       [[0.44633162, 0.60684282, 0.96393795],
        [0.13904548, 0.73709456, 0.37711571],
        [0.79070871, 0.15560342, 0.12328662],
        [0.40389155, 0.97863405, 0.0114599 ],
        [0.05951757, 0.10889081, 0.85272281],
        [0.71514612, 0.52312704, 0.59877454],
        [0.64931396, 0.23505158, 0.48710589],
        [0.30097968, 0.36957094, 0.23358374],
        [0.99369356, 0.78380717, 0.74090153]]])
>>> Faure(3,randomize="LMS DS",seed=7,replications=2)(9)
array([[[1.90132648e-01, 7.86261604e-01, 1.43261685e-01],
        [4.27494577e-01, 6.21249393e-01, 4.43831847e-01],
        [9.93922386e-01, 9.24890107e-02, 9.54261973e-01],
        [8.29146531e-02, 4.49180940e-01, 7.37185290e-01],
        [6.53610094e-01, 2.95352182e-01, 2.47631786e-01],
        [8.86856776e-01, 7.55465004e-01, 5.48200701e-01],
        [2.96893357e-01, 1.12311024e-01, 6.56285230e-01],
        [5.34102847e-01, 9.58324412e-01, 8.45725185e-01],
        [7.67349545e-01, 4.29366454e-01, 3.51846617e-02]],

       [[9.98454735e-01, 2.53648344e-01, 7.13434357e-01],
        [3.05609133e-01, 4.45718366e-01, 6.36691290e-01],
        [6.54114646e-01, 6.91349078e-01, 1.11471149e-01],
        [7.11691963e-01, 6.25495629e-01, 2.23136805e-01],
        [5.60878749e-02, 8.17408217e-01, 8.25559119e-01],
        [4.00275568e-01, 5.89230529e-02, 4.15633066e-01],
        [7.91236327e-01, 9.93070231e-01, 5.25932625e-01],
        [1.39697286e-01, 1.89249857e-01, 7.14552095e-04],
        [4.46998476e-01, 4.30617921e-01, 9.36315716e-01]]])
>>> Faure(3,randomize="LMS",seed=7,replications=2)(9,warn=False)
array([[[0.        , 0.        , 0.        ],
        [0.5706961 , 0.84756013, 0.63390413],
        [0.80394278, 0.65243798, 0.82336297],
        [0.22678173, 0.67547662, 0.93144564],
        [0.46399184, 0.50916041, 0.12088811],
        [0.69723852, 0.31536484, 0.40911072],
        [0.11956371, 0.33845241, 0.52953528],
        [0.35692624, 0.18443172, 0.70664719],
        [0.92335406, 0.97711587, 0.22942236]],

       [[0.        , 0.        , 0.        ],
        [0.34851179, 0.57896969, 0.93575568],
        [0.69285188, 0.80808804, 0.52582983],
        [0.20726949, 0.37184791, 0.6375009 ],
        [0.55146283, 0.94654367, 0.11228285],
        [0.8628786 , 0.17977726, 0.71455089],
        [0.2908789 , 0.74353714, 0.82500264],
        [0.59818511, 0.31838531, 0.41492625],
        [0.94242359, 0.54735751, 0.22303965]]])
>>> Faure(3,randomize="DS",seed=7,replications=2)(9)
array([[[0.68058188, 0.90130449, 0.36582045],
        [0.01391521, 0.23463782, 0.69915378],
        [0.34724854, 0.56797116, 0.03248711],
        [0.79169299, 0.0124156 , 0.14359822],
        [0.12502632, 0.34574893, 0.47693156],
        [0.45835966, 0.67908227, 0.81026489],
        [0.9028041 , 0.45686005, 0.921376  ],
        [0.23613743, 0.79019338, 0.25470933],
        [0.56947077, 0.12352671, 0.58804267]],

       [[0.95466555, 0.56062266, 0.29687589],
        [0.28799888, 0.89395599, 0.63020922],
        [0.62133222, 0.22728933, 0.96354256],
        [0.73244333, 0.67173377, 0.74132033],
        [0.06577666, 0.0050671 , 0.07465367],
        [0.39911   , 0.33840044, 0.407987  ],
        [0.84355444, 0.11617822, 0.51909811],
        [0.17688777, 0.44951155, 0.85243145],
        [0.51022111, 0.78284488, 0.18576478]]])
>>> Faure(3,randomize="DP",seed=7,replications=2)(9)
array([[[0.62331102, 0.65595286, 0.5208953 ],
        [0.28997769, 0.98928619, 0.18756197],
        [0.95664436, 0.32261952, 0.85422863],
        [0.51219991, 0.87817508, 0.96533975],
        [0.17886658, 0.21150841, 0.63200641],
        [0.84553325, 0.54484175, 0.29867308],
        [0.4010888 , 0.1003973 , 0.07645086],
        [0.06775547, 0.43373064, 0.74311752],
        [0.73442213, 0.76706397, 0.40978419]],

       [[0.65992302, 0.22190291, 0.57359374],
        [0.32658969, 0.88856957, 0.90692707],
        [0.99325635, 0.55523624, 0.2402604 ],
        [0.54881191, 0.77745846, 0.12914929],
        [0.21547857, 0.44412513, 0.46248263],
        [0.88214524, 0.1107918 , 0.79581596],
        [0.4377008 , 0.66634735, 0.68470485],
        [0.10436746, 0.33301402, 0.01803818],
        [0.77103413, 0.99968068, 0.35137152]]])
>>> Faure(3,randomize="NUS",seed=7,replications=2)(9)
array([[[0.077534  , 0.99285569, 0.64186774],
        [0.82225806, 0.4963679 , 0.86453228],
        [0.48498632, 0.0833722 , 0.06011164],
        [0.15691802, 0.60236449, 0.11226224],
        [0.95445636, 0.2796012 , 0.54752025],
        [0.61890325, 0.84183901, 0.70654863],
        [0.25089638, 0.17805972, 0.95988146],
        [0.68344029, 0.77065782, 0.26676153],
        [0.4322891 , 0.40799837, 0.34911626]],

       [[0.03739061, 0.16158904, 0.85379075],
        [0.7635801 , 0.40075819, 0.56644152],
        [0.49937232, 0.74097216, 0.02663559],
        [0.24837929, 0.59622709, 0.12164937],
        [0.79067346, 0.83319946, 0.67040549],
        [0.40808716, 0.31605175, 0.51589378],
        [0.15774296, 0.96060354, 0.35917763],
        [0.96867006, 0.05032588, 0.24759375],
        [0.59326363, 0.50120469, 0.9906825 ]]])

Higher order Faure

>>> Faure(3,randomize="LMS DP",seed=7,alpha=2)(9)
array([[0.07060326, 0.24965078, 0.49971375],
       [0.9104272 , 0.77359118, 0.02813304],
       [0.52709818, 0.51562414, 0.97220901],
       [0.23236901, 0.55967785, 0.28909699],
       [0.88746566, 0.08382775, 0.77956468],
       [0.39256201, 0.81762837, 0.43144077],
       [0.20940728, 0.91289379, 0.71118926],
       [0.71589223, 0.4203758 , 0.56884549],
       [0.59136462, 0.16673034, 0.219807  ]])
>>> Faure(3,randomize="LMS DS",seed=7,alpha=2)(9)
array([[0.27004662, 0.64652898, 0.04873353],
       [0.35569696, 0.06066804, 0.54773651],
       [0.88670561, 0.79188841, 0.90363807],
       [0.03415605, 0.99796556, 0.62789873],
       [0.5267559 , 0.39579989, 0.76064852],
       [0.94711696, 0.14327166, 0.11150335],
       [0.20948284, 0.29996665, 0.83571342],
       [0.62847216, 0.71019873, 0.30272608],
       [0.67845718, 0.45371208, 0.36140179]])
>>> Faure(3,randomize="LMS",seed=7,alpha=2)(9,warn=False)
array([[0.        , 0.        , 0.        ],
       [0.53148918, 0.85880934, 0.51685452],
       [0.95184459, 0.60232457, 0.98293576],
       [0.21282738, 0.35205188, 0.58389855],
       [0.59616274, 0.20653524, 0.72441959],
       [0.68272695, 0.94187019, 0.19189138],
       [0.2733855 , 0.7034865 , 0.79264466],
       [0.35996672, 0.5457835 , 0.27107157],
       [0.85440183, 0.28913878, 0.43628397]])
>>> Faure(3,randomize="NUS",seed=7,alpha=2)(9)
array([[0.25822173, 0.41712395, 0.27892432],
       [0.88275133, 0.97811616, 0.46301445],
       [0.3846734 , 0.21568889, 0.71276888],
       [0.16342299, 0.03041244, 0.40224624],
       [0.71060873, 0.62350799, 0.97150042],
       [0.57999438, 0.78144529, 0.21321557],
       [0.10439042, 0.71635743, 0.78194125],
       [0.89714231, 0.26123842, 0.01970874],
       [0.52510616, 0.46613858, 0.65895473]])
>>> Faure(3,randomize="False",seed=7,alpha=2)(9,warn=False)
array([[0.        , 0.        , 0.        ],
       [0.44444444, 0.44444444, 0.44444444],
       [0.88888889, 0.88888889, 0.88888889],
       [0.16049383, 0.71604938, 0.60493827],
       [0.60493827, 0.16049383, 0.71604938],
       [0.71604938, 0.60493827, 0.16049383],
       [0.32098765, 0.43209877, 0.87654321],
       [0.43209877, 0.87654321, 0.32098765],
       [0.87654321, 0.32098765, 0.43209877]])

Replications of higher order Faure

>>> Faure(3,randomize="LMS DP",seed=7,alpha=2,replications=2)(9)
array([[[0.65006542, 0.84004771, 0.39377772],
        [0.73541117, 0.25289783, 0.11639162],
        [0.11843564, 0.44363425, 0.98988211],
        [0.37561943, 0.50110654, 0.26965135],
        [0.79429224, 0.90984833, 0.69823833],
        [0.32592209, 0.08813056, 0.53196334],
        [0.47565296, 0.15886247, 0.84891661],
        [0.96885711, 0.55947605, 0.57425894],
        [0.05559154, 0.74599626, 0.07691998]],

       [[0.01647589, 0.79799093, 0.98866597],
        [0.79175933, 0.39051971, 0.1440758 ],
        [0.57923112, 0.31651741, 0.39977363],
        [0.19716729, 0.14023077, 0.0888261 ],
        [0.9681839 , 0.73296284, 0.5775844 ],
        [0.41014543, 0.6588137 , 0.83313225],
        [0.28590138, 0.45066659, 0.52218222],
        [0.74006249, 0.04320102, 0.67774517],
        [0.51107317, 0.96909703, 0.26664273]]])
>>> Faure(3,randomize="LMS DS",seed=7,alpha=2,replications=2)(9)
array([[[0.82257501, 0.14639929, 0.60116121],
        [0.24249012, 0.98867933, 0.77186064],
        [0.43925164, 0.40287298, 0.12683077],
        [0.91956054, 0.44754172, 0.85119388],
        [0.08067361, 0.30216242, 0.31317864],
        [0.49966237, 0.71234433, 0.33572993],
        [0.75789644, 0.7949646 , 0.05999059],
        [0.17827381, 0.65360268, 0.52607183],
        [0.55976324, 0.05143265, 0.91398251]],

       [[0.90790039, 0.30655016, 0.02228959],
        [0.34969458, 0.71397805, 0.9779997 ],
        [0.13261487, 0.45477877, 0.50021876],
        [0.71256557, 0.66536938, 0.789166  ],
        [0.48786265, 0.07285372, 0.41137525],
        [0.26253474, 0.81344564, 0.2669434 ],
        [0.87953962, 0.98487119, 0.58886448],
        [0.6628982 , 0.39219753, 0.2112249 ],
        [0.10438938, 0.1329926 , 0.73328965]]])
>>> Faure(3,randomize="LMS",seed=7,alpha=2,replications=2)(9,warn=False)
array([[[0.        , 0.        , 0.        ],
        [0.53148918, 0.85880934, 0.51685452],
        [0.95184459, 0.60232457, 0.98293576],
        [0.21282738, 0.35205188, 0.58389855],
        [0.59616274, 0.20653524, 0.72441959],
        [0.68272695, 0.94187019, 0.19189138],
        [0.2733855 , 0.7034865 , 0.79264466],
        [0.35996672, 0.5457835 , 0.27107157],
        [0.85440183, 0.28913878, 0.43628397]],

       [[0.        , 0.        , 0.        ],
        [0.78765092, 0.85192948, 0.96805852],
        [0.56277838, 0.59274149, 0.49028114],
        [0.18069211, 0.36070495, 0.80396678],
        [0.95171191, 0.21257791, 0.43852003],
        [0.40600476, 0.9532375 , 0.29409298],
        [0.31884304, 0.68248734, 0.57074159],
        [0.76064476, 0.53425869, 0.20544726],
        [0.53167411, 0.27502561, 0.72751996]]])
>>> Faure(3,randomize="NUS",seed=7,alpha=2,replications=2)(9)
array([[[0.25822173, 0.41712395, 0.27892432],
        [0.88275133, 0.97811616, 0.46301445],
        [0.3846734 , 0.21568889, 0.71276888],
        [0.16342299, 0.03041244, 0.40224624],
        [0.71060873, 0.62350799, 0.97150042],
        [0.57999438, 0.78144529, 0.21321557],
        [0.10439042, 0.71635743, 0.78194125],
        [0.89714231, 0.26123842, 0.01970874],
        [0.52510616, 0.46613858, 0.65895473]],

       [[0.1950677 , 0.52270739, 0.50736994],
        [0.95117968, 0.09639097, 0.23235875],
        [0.37944055, 0.98868015, 0.6693519 ],
        [0.27700817, 0.83696832, 0.072855  ],
        [0.6981059 , 0.3763676 , 0.8630192 ],
        [0.46377485, 0.28179682, 0.64253665],
        [0.01958311, 0.13843641, 0.9588172 ],
        [0.87752263, 0.7029967 , 0.35134227],
        [0.64185819, 0.55907117, 0.19929854]]])
Source code in qmcpy/discrete_distribution/digital_net_any_bases/digital_net_any_bases.py 
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
def __init__(self,
             dimension = 1,
             replications = None,
             seed = None, 
             randomize = 'LMS DP',
             bases_generating_matrices = None,
             t = None,
             alpha = 1,
             n_lim = 2**32,
             warn = True):
    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. 
        bases_generating_matrices (Union[str, tuple]: Specify the bases and the generating matrices.

            - `"HALTON"` will use Halton generating matrices.
            - `"FAURE"` will use Faure generating matrices .
            - `bases,generating_matrices` requires 

                - `bases` is an `np.ndarray` of integers with shape $(,d)$ or $(r,d)$ where $d$ is the number of dimensions and $r$ is the number of replications.
                - `generating_matrices` is an `np.ndarray` of integers with shape $(d,m_\mathrm{max},t_\mathrm{max})$ or $(r,d,m_\mathrm{max},t_\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.

        t (int): Number of digits *after* randomization. The number of digits in the generating matrices is inferred.
        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.  
        n_lim (int): Maximum number of compatible points, determines the number of rows in the generating matrices. 
        warn (bool): If `False`, suppress warnings in construction 
    """
    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)
    if bases_generating_matrices is None:
        if self.DEFAULT_GENERATING_MATRICES=="HALTON":
            self.type_bases_generating_matrices = "HALTON"
            d_limit = len(self.all_primes)
        elif self.DEFAULT_GENERATING_MATRICES=="FAURE":
            self.type_bases_generating_matrices = "FAURE"
            d_limit = int(self.all_primes[-1])
        else:
            raise ParameterError("must supply bases_generating_matrices")
    else:
        self.type_bases_generating_matrices = "CUSTOM"
        assert len(bases_generating_matrices)==2
        bases,generating_matrices = bases_generating_matrices
        assert isinstance(generating_matrices,np.ndarray)
        assert generating_matrices.ndim==3 or generating_matrices.ndim==4 
        d_limit = generating_matrices.shape[1]
        if np.isscalar(bases):
            assert bases>0
            assert bases%1==0 
            bases = int(bases)*np.ones(d_limit,dtype=int)
        assert bases.ndim==1 or bases.ndim==2
    self.input_t = deepcopy(t) 
    self.input_bases_generating_matrices = deepcopy(bases_generating_matrices)
    super(DigitalNetAnyBases,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":
        assert self.type_bases_generating_matrices=="HALTON", "QRNG randomization is only applicable for the Halton generator."
        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.alpha = alpha
    assert self.alpha>=1
    assert self.alpha%1==0
    if self.alpha>1:
        assert (self.dvec==np.arange(self.d)).all(), "digital interlacing requires dimension is an int"
    self.dtalpha = self.alpha*self.d 
    if self.type_bases_generating_matrices=="HALTON":
        self.bases = self.all_primes[self.dvec][None,:]
        self.m_max = int(np.ceil(np.log(self.n_limit)/np.log(self.bases.min())))
        self._t_curr = self.m_max
        if t is None: t = int(np.ceil(-np.log(2**(-63))/np.log(self.bases.min())))
        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))
    elif self.type_bases_generating_matrices=="FAURE":
        assert (self.dvec==np.arange(self.d)).all(), "Faure requires dimension is an int"
        p = self.all_primes[np.argmax(self.all_primes>=self.d)]
        self.bases = p*np.ones((1,self.dtalpha),dtype=np.uint64)
        self.m_max = int(np.ceil(np.log(self.n_limit)/np.log(p)))
        if t is None: t = int(np.ceil(-np.log(2**(-63))/np.log(p)))
        self._t_curr = self.m_max
        self.t = self.m_max if self.m_max>t else t
        self.C = np.ones((self.dtalpha,1,1),dtype=np.uint64)*np.eye(self._t_curr,dtype=np.uint64)
        if self.dtalpha>1:
            for a in range(self._t_curr):
                for b in range(a+1):
                    self.C[1,a,b] = comb(a,b)%p
        if self.dtalpha>2:
            for k in range(2,self.dtalpha):
                for a in range(self._t_curr):
                    for b in range(a+1):
                        self.C[k,a,b] = (int(self.C[1,a,b])*((k**(a-b))%p))%p
        self.C = self.C[None,:,:,:]
        self.faure_base = p
    else:
        self.bases = bases.astype(np.uint64)
        if self.bases.ndim==1: self.bases = self.bases[None,:]
        assert self.bases.shape[1]>=self.dtalpha
        if self.alpha==1:
            self.bases = self.bases[:,self.dvec]
        else:
            self.bases = self.bases[:,:self.dtalpha]
        self.C = generating_matrices.astype(np.uint64)
        if self.C.ndim==3: self.C = self.C[None,:,:,:]
        assert self.C.shape[1]>=self.dtalpha
        if self.alpha==1:
            self.C = self.C[:,self.dvec,:,:]
        else:
            self.C = self.C[:,:self.dtalpha,:,:]
        self.m_max,self._t_curr = self.C.shape[-2:]
        if t is None: t = int(np.ceil(-np.log(2**(-63))/np.log(self.bases.min())))
        self.t = self.m_max if self.m_max>t else t
        assert (0<=self.C).all()
        assert (self.C<self.bases[:,:,None,None]).all()
        self.C = np.ascontiguousarray(self.C) 
        self.bases = np.ascontiguousarray(self.bases)
    if self.alpha>1:
        assert (self.bases==self.bases[0,0]).all(), "alpha>1 performs digital interlacing which requires the same base across dimensions and replications."
        if warn and self.m_max!=self._t_curr:
            warnings.warn("Digital interlacing is often performed on generating matrices with the number of columns (m_max = %d) equal to the number of rows (_t_curr = %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)
    assert self.bases.ndim==2
    assert self.bases.shape[-1]==self.dtalpha
    assert self.bases.shape[0]==1 or self.bases.shape[0]==self.replications
    assert self.C.ndim==4
    assert self.C.shape[-3:]==(self.dtalpha,self.m_max,self._t_curr)
    assert self.C.shape[0]==1 or self.C.shape[0]==self.replications
    r_b = self.bases.shape[0]
    r_C = self.C.shape[0]
    if self.randomize=="FALSE":
        if self.alpha>1:
            t_alpha = min(self.t,self._t_curr*self.alpha)
            C_ho = np.empty((self.replications,self.d,self.m_max,t_alpha),dtype=np.uint64)
            qmctoolscl.gdn_interlace(np.uint64(self.replications),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),self.C,C_ho)
            self.C = C_ho
            self._t_curr = t_alpha
            self.t = t_alpha
            self.bases = self.bases[:,:self.d]
    elif self.randomize=="DP":
        self.perms = qmctoolscl.gdn_get_digital_permutations(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(r_b),self.bases)
    elif self.randomize=="DS":
        self.rshift = qmctoolscl.gdn_get_digital_shifts(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(r_b),self.bases)
    elif self.randomize in ["LMS","LMS DS","LMS DP"]:
        if self.alpha==1:
            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(r_b),self.bases)
            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(r_C),np.uint64(r_b),np.uint64(self._t_curr),np.uint64(self.t),self.bases,S,self.C,C_lms,backend="c")
            self.C = C_lms
            self._t_curr = self.t
        else:
            t_dig = np.ceil(max(self.t/self.alpha,self._t_curr))
            S = qmctoolscl.gdn_get_linear_scramble_matrix(self.rng,np.uint64(self.replications),np.uint64(self.dtalpha),np.uint64(self._t_curr),np.uint64(t_dig),np.uint64(r_b),self.bases)
            C_lms = np.empty((self.replications,self.dtalpha,self.m_max,self.t),dtype=np.uint64)
            qmctoolscl.gdn_linear_matrix_scramble(np.uint64(self.replications),np.uint64(self.dtalpha),np.uint64(self.m_max),np.uint64(r_C),np.uint64(r_b),np.uint64(self._t_curr),np.uint64(t_dig),self.bases,S,self.C,C_lms,backend="c")
            C_lms_ho = np.empty((self.replications,self.d,self.m_max,self.t),dtype=np.uint64)
            qmctoolscl.gdn_interlace(np.uint64(self.replications),np.uint64(self.d),np.uint64(self.m_max),np.uint64(self.dtalpha),np.uint64(t_dig),np.uint64(self.t),np.uint64(self.alpha),C_lms,C_lms_ho)
            self.C = C_lms_ho
            self._t_curr = self.t
            self.t = self._t_curr
            self.bases = self.bases[:,:self.d]
        if self.randomize=="LMS DP":
            self.perms = qmctoolscl.gdn_get_digital_permutations(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(r_b),self.bases)
        elif self.randomize=="LMS DS":
            self.rshift = qmctoolscl.gdn_get_digital_shifts(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(r_b),self.bases)
    elif self.randomize=="NUS":
        if self.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_gdn() for i in range(self.replications*self.d)]).reshape(self.replications,self.d)
        else:
            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_gdn() for i in range(self.replications*self.dtalpha)]).reshape(self.replications,self.dtalpha)
    assert self.C.ndim==4 and (self.C.shape[0]==1 or self.C.shape[0]==self.replications) and self.C.shape[1]==(self.dtalpha if self.randomize=="NUS" else self.d) and self.C.shape[2]==self.m_max and self.C.shape[3]==self._t_curr
    assert self.bases.ndim==2 and (self.bases.shape[0]==1 or self.bases.shape[0]==self.replications) and self.bases.shape[1]==(self.dtalpha if self.randomize=="NUS" else self.d)
    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])
    if warn and (self.bases==2).all():
        warnings.warn("It is more efficient to use DigitalNetB2 instead of DigitalNetAnyBases when all bases are 2")
    self.warn = warn

DigitalNetAnyBases

Bases: AbstractLDDiscreteDistribution

Low discrepancy digital net with arbitrary bases for each dimension.

Note
  • Digital net samples sizes should be products of powers of bases, i.e., a digital net with bases \((b_1,\dots,b_d)\) will prefer sample sizes \(n = b_1^{p_1} \cdots b_d^{p_d}\) for some \(p_1,\dots,p_d \in \mathbb{N}_0\).
  • The first point of an unrandomized digital net is the origin.

  • The construction of higher order digital nets requires the same base for each dimension. To construct 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.

A few examples below showcase how to pass in custom bases and generating matrices. Many other examples can be found in the Halton and Faure implementations

Examples:

>>> bases = 3 
>>> generating_matrices = np.array(
...     [
...         [[1, 0, 0],
...          [0, 1, 0],
...          [1, 2, 1]],
...         [[2, 0, 0],
...          [1, 2, 0],
...          [0, 2, 2]]
...     ],
...     dtype=np.uint64)
>>> DigitalNetAnyBases(2,randomize="False",bases_generating_matrices=(bases,generating_matrices))(9,warn=False)
array([[0.        , 0.        ],
       [0.33333333, 0.66666667],
       [0.66666667, 0.33333333],
       [0.11111111, 0.55555556],
       [0.44444444, 0.22222222],
       [0.77777778, 0.88888889],
       [0.22222222, 0.77777778],
       [0.55555556, 0.44444444],
       [0.88888889, 0.11111111]])
>>> DigitalNetAnyBases(1,randomize="False",alpha=2,bases_generating_matrices=(bases,generating_matrices))(9,warn=False)
array([[0.        ],
       [0.55555556],
       [0.77777778],
       [0.17283951],
       [0.39506173],
       [0.95061728],
       [0.30864198],
       [0.5308642 ],
       [0.75308642]])

>>> rng = np.random.Generator(np.random.PCG64(7))
>>> bases = np.array(
...     [[2,5,7,23],
...      [3,11,13,17]],
...     dtype=np.uint64)
>>> generating_matrices = np.stack(
...     [
...         np.stack(
...             [   np.tril(rng.integers(0,2,(10,10))),
...                 np.tril(rng.integers(0,5,(10,10))),
...                 np.tril(rng.integers(0,7,(10,10))),
...                 np.tril(rng.integers(0,23,(10,10))),],axis=0),
...         np.stack(
...             [   np.tril(rng.integers(0,3,(10,10))),
...                 np.tril(rng.integers(0,11,(10,10))),
...                 np.tril(rng.integers(0,13,(10,10))),
...                 np.tril(rng.integers(0,17,(10,10))),],axis=0),
...     ],axis=0).astype(np.uint64)
>>> discrete_distrib = DigitalNetAnyBases(
...     dimension = 4,
...     randomize = "NUS",
...     seed = 7,
...     replications = 2,
...     bases_generating_matrices=(bases,generating_matrices),
...     )
>>> discrete_distrib(10,warn=False)
array([[[0.141964  , 0.50711584, 0.75766621, 0.43359373],
        [0.65536579, 0.66327673, 0.21825902, 0.93901529],
        [0.46955206, 0.21554445, 0.54265048, 0.17421341],
        [0.78505432, 0.99432677, 0.08408121, 0.8317411 ],
        [0.65536579, 0.05422203, 0.37171211, 0.05616883],
        [0.141964  , 0.21554445, 0.94198662, 0.89162246],
        [0.78505432, 0.99432677, 0.62395426, 0.70351649],
        [0.46955206, 0.05422203, 0.51786701, 0.75215221],
        [0.98479561, 0.50711584, 0.06267682, 0.59165814],
        [0.31524867, 0.66327673, 0.29204202, 0.29200658]],

       [[0.16158904, 0.00911626, 0.99631587, 0.7984386 ],
        [0.33578478, 0.98734618, 0.48165571, 0.84189861],
        [0.84596814, 0.67093277, 0.72532656, 0.49339621],
        [0.33578478, 0.54281197, 0.26069214, 0.43832903],
        [0.84596814, 0.58529648, 0.80521628, 0.15910404],
        [0.16158904, 0.4486066 , 0.84973668, 0.6750655 ],
        [0.84596814, 0.21294277, 0.40387763, 0.19059409],
        [0.16158904, 0.74257456, 0.19604142, 0.98366484],
        [0.33578478, 0.31746094, 0.35948446, 0.75911922],
        [0.64593022, 0.11007165, 0.63174328, 0.55910368]]])

>>> bases = 2
>>> generating_matrices = np.array([
...     [[1, 0, 0],
...      [0, 1, 0],
...      [1, 0, 1]],
...     [[1, 0, 0],
...      [1, 1, 0],
...      [0, 1, 1]]],dtype=np.uint64)
>>> x = DigitalNetAnyBases(2,randomize="False",bases_generating_matrices=(bases,generating_matrices),warn=False)(8,warn=False)
>>> x
array([[0.   , 0.   ],
       [0.5  , 0.5  ],
       [0.25 , 0.75 ],
       [0.75 , 0.25 ],
       [0.625, 0.375],
       [0.125, 0.875],
       [0.875, 0.625],
       [0.375, 0.125]])
>>> generating_matrices_b2 = (generating_matrices*2**np.array([2,1,0])).sum(-1).astype(np.uint64)
>>> x_b2 = DigitalNetB2(2,randomize="False",order="radical inverse",generating_matrices=generating_matrices_b2,t=3,msb=True)(8,warn=False)
>>> x_b2
array([[0.   , 0.   ],
       [0.5  , 0.5  ],
       [0.25 , 0.75 ],
       [0.75 , 0.25 ],
       [0.625, 0.375],
       [0.125, 0.875],
       [0.875, 0.625],
       [0.375, 0.125]])
>>> bool((x==x_b2).all())
True

References:

  1. Dick, Josef, and Friedrich Pillichshammer.
    Digital nets and sequences: discrepancy theory and quasi–Monte Carlo integration.
    Cambridge University Press, 2010.

  2. Sorokin, Aleksei.
    "QMCPy: A Python Software for Randomized Low-Discrepancy Sequences, Quasi-Monte Carlo, and Fast Kernel Methods"
    arXiv preprint arXiv:2502.14256 (2025).

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'
bases_generating_matrices (Union[str, tuple]

Specify the bases and the generating matrices.

  • "HALTON" will use Halton generating matrices.
  • "FAURE" will use Faure generating matrices .

  • bases,generating_matrices requires

    • bases is an np.ndarray of integers with shape \((,d)\) or \((r,d)\) where \(d\) is the number of dimensions and \(r\) is the number of replications.
    • generating_matrices is an np.ndarray of integers with shape \((d,m_\mathrm{max},t_\mathrm{max})\) or \((r,d,m_\mathrm{max},t_\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.
 required
t int

Number of digits after randomization. The number of digits in the generating matrices is inferred.

 None
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
n_lim int

Maximum number of compatible points, determines the number of rows in the generating matrices.

2 ** 32
warn bool

If False, suppress warnings in construction

 True
Source code in qmcpy/discrete_distribution/digital_net_any_bases/digital_net_any_bases.py 
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
def __init__(self,
             dimension = 1,
             replications = None,
             seed = None, 
             randomize = 'LMS DP',
             bases_generating_matrices = None,
             t = None,
             alpha = 1,
             n_lim = 2**32,
             warn = True):
    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. 
        bases_generating_matrices (Union[str, tuple]: Specify the bases and the generating matrices.

            - `"HALTON"` will use Halton generating matrices.
            - `"FAURE"` will use Faure generating matrices .
            - `bases,generating_matrices` requires 

                - `bases` is an `np.ndarray` of integers with shape $(,d)$ or $(r,d)$ where $d$ is the number of dimensions and $r$ is the number of replications.
                - `generating_matrices` is an `np.ndarray` of integers with shape $(d,m_\mathrm{max},t_\mathrm{max})$ or $(r,d,m_\mathrm{max},t_\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.

        t (int): Number of digits *after* randomization. The number of digits in the generating matrices is inferred.
        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.  
        n_lim (int): Maximum number of compatible points, determines the number of rows in the generating matrices. 
        warn (bool): If `False`, suppress warnings in construction 
    """
    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)
    if bases_generating_matrices is None:
        if self.DEFAULT_GENERATING_MATRICES=="HALTON":
            self.type_bases_generating_matrices = "HALTON"
            d_limit = len(self.all_primes)
        elif self.DEFAULT_GENERATING_MATRICES=="FAURE":
            self.type_bases_generating_matrices = "FAURE"
            d_limit = int(self.all_primes[-1])
        else:
            raise ParameterError("must supply bases_generating_matrices")
    else:
        self.type_bases_generating_matrices = "CUSTOM"
        assert len(bases_generating_matrices)==2
        bases,generating_matrices = bases_generating_matrices
        assert isinstance(generating_matrices,np.ndarray)
        assert generating_matrices.ndim==3 or generating_matrices.ndim==4 
        d_limit = generating_matrices.shape[1]
        if np.isscalar(bases):
            assert bases>0
            assert bases%1==0 
            bases = int(bases)*np.ones(d_limit,dtype=int)
        assert bases.ndim==1 or bases.ndim==2
    self.input_t = deepcopy(t) 
    self.input_bases_generating_matrices = deepcopy(bases_generating_matrices)
    super(DigitalNetAnyBases,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":
        assert self.type_bases_generating_matrices=="HALTON", "QRNG randomization is only applicable for the Halton generator."
        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.alpha = alpha
    assert self.alpha>=1
    assert self.alpha%1==0
    if self.alpha>1:
        assert (self.dvec==np.arange(self.d)).all(), "digital interlacing requires dimension is an int"
    self.dtalpha = self.alpha*self.d 
    if self.type_bases_generating_matrices=="HALTON":
        self.bases = self.all_primes[self.dvec][None,:]
        self.m_max = int(np.ceil(np.log(self.n_limit)/np.log(self.bases.min())))
        self._t_curr = self.m_max
        if t is None: t = int(np.ceil(-np.log(2**(-63))/np.log(self.bases.min())))
        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))
    elif self.type_bases_generating_matrices=="FAURE":
        assert (self.dvec==np.arange(self.d)).all(), "Faure requires dimension is an int"
        p = self.all_primes[np.argmax(self.all_primes>=self.d)]
        self.bases = p*np.ones((1,self.dtalpha),dtype=np.uint64)
        self.m_max = int(np.ceil(np.log(self.n_limit)/np.log(p)))
        if t is None: t = int(np.ceil(-np.log(2**(-63))/np.log(p)))
        self._t_curr = self.m_max
        self.t = self.m_max if self.m_max>t else t
        self.C = np.ones((self.dtalpha,1,1),dtype=np.uint64)*np.eye(self._t_curr,dtype=np.uint64)
        if self.dtalpha>1:
            for a in range(self._t_curr):
                for b in range(a+1):
                    self.C[1,a,b] = comb(a,b)%p
        if self.dtalpha>2:
            for k in range(2,self.dtalpha):
                for a in range(self._t_curr):
                    for b in range(a+1):
                        self.C[k,a,b] = (int(self.C[1,a,b])*((k**(a-b))%p))%p
        self.C = self.C[None,:,:,:]
        self.faure_base = p
    else:
        self.bases = bases.astype(np.uint64)
        if self.bases.ndim==1: self.bases = self.bases[None,:]
        assert self.bases.shape[1]>=self.dtalpha
        if self.alpha==1:
            self.bases = self.bases[:,self.dvec]
        else:
            self.bases = self.bases[:,:self.dtalpha]
        self.C = generating_matrices.astype(np.uint64)
        if self.C.ndim==3: self.C = self.C[None,:,:,:]
        assert self.C.shape[1]>=self.dtalpha
        if self.alpha==1:
            self.C = self.C[:,self.dvec,:,:]
        else:
            self.C = self.C[:,:self.dtalpha,:,:]
        self.m_max,self._t_curr = self.C.shape[-2:]
        if t is None: t = int(np.ceil(-np.log(2**(-63))/np.log(self.bases.min())))
        self.t = self.m_max if self.m_max>t else t
        assert (0<=self.C).all()
        assert (self.C<self.bases[:,:,None,None]).all()
        self.C = np.ascontiguousarray(self.C) 
        self.bases = np.ascontiguousarray(self.bases)
    if self.alpha>1:
        assert (self.bases==self.bases[0,0]).all(), "alpha>1 performs digital interlacing which requires the same base across dimensions and replications."
        if warn and self.m_max!=self._t_curr:
            warnings.warn("Digital interlacing is often performed on generating matrices with the number of columns (m_max = %d) equal to the number of rows (_t_curr = %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)
    assert self.bases.ndim==2
    assert self.bases.shape[-1]==self.dtalpha
    assert self.bases.shape[0]==1 or self.bases.shape[0]==self.replications
    assert self.C.ndim==4
    assert self.C.shape[-3:]==(self.dtalpha,self.m_max,self._t_curr)
    assert self.C.shape[0]==1 or self.C.shape[0]==self.replications
    r_b = self.bases.shape[0]
    r_C = self.C.shape[0]
    if self.randomize=="FALSE":
        if self.alpha>1:
            t_alpha = min(self.t,self._t_curr*self.alpha)
            C_ho = np.empty((self.replications,self.d,self.m_max,t_alpha),dtype=np.uint64)
            qmctoolscl.gdn_interlace(np.uint64(self.replications),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),self.C,C_ho)
            self.C = C_ho
            self._t_curr = t_alpha
            self.t = t_alpha
            self.bases = self.bases[:,:self.d]
    elif self.randomize=="DP":
        self.perms = qmctoolscl.gdn_get_digital_permutations(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(r_b),self.bases)
    elif self.randomize=="DS":
        self.rshift = qmctoolscl.gdn_get_digital_shifts(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(r_b),self.bases)
    elif self.randomize in ["LMS","LMS DS","LMS DP"]:
        if self.alpha==1:
            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(r_b),self.bases)
            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(r_C),np.uint64(r_b),np.uint64(self._t_curr),np.uint64(self.t),self.bases,S,self.C,C_lms,backend="c")
            self.C = C_lms
            self._t_curr = self.t
        else:
            t_dig = np.ceil(max(self.t/self.alpha,self._t_curr))
            S = qmctoolscl.gdn_get_linear_scramble_matrix(self.rng,np.uint64(self.replications),np.uint64(self.dtalpha),np.uint64(self._t_curr),np.uint64(t_dig),np.uint64(r_b),self.bases)
            C_lms = np.empty((self.replications,self.dtalpha,self.m_max,self.t),dtype=np.uint64)
            qmctoolscl.gdn_linear_matrix_scramble(np.uint64(self.replications),np.uint64(self.dtalpha),np.uint64(self.m_max),np.uint64(r_C),np.uint64(r_b),np.uint64(self._t_curr),np.uint64(t_dig),self.bases,S,self.C,C_lms,backend="c")
            C_lms_ho = np.empty((self.replications,self.d,self.m_max,self.t),dtype=np.uint64)
            qmctoolscl.gdn_interlace(np.uint64(self.replications),np.uint64(self.d),np.uint64(self.m_max),np.uint64(self.dtalpha),np.uint64(t_dig),np.uint64(self.t),np.uint64(self.alpha),C_lms,C_lms_ho)
            self.C = C_lms_ho
            self._t_curr = self.t
            self.t = self._t_curr
            self.bases = self.bases[:,:self.d]
        if self.randomize=="LMS DP":
            self.perms = qmctoolscl.gdn_get_digital_permutations(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(r_b),self.bases)
        elif self.randomize=="LMS DS":
            self.rshift = qmctoolscl.gdn_get_digital_shifts(self.rng,np.uint64(self.replications),np.uint64(self.d),self.t,np.uint64(r_b),self.bases)
    elif self.randomize=="NUS":
        if self.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_gdn() for i in range(self.replications*self.d)]).reshape(self.replications,self.d)
        else:
            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_gdn() for i in range(self.replications*self.dtalpha)]).reshape(self.replications,self.dtalpha)
    assert self.C.ndim==4 and (self.C.shape[0]==1 or self.C.shape[0]==self.replications) and self.C.shape[1]==(self.dtalpha if self.randomize=="NUS" else self.d) and self.C.shape[2]==self.m_max and self.C.shape[3]==self._t_curr
    assert self.bases.ndim==2 and (self.bases.shape[0]==1 or self.bases.shape[0]==self.replications) and self.bases.shape[1]==(self.dtalpha if self.randomize=="NUS" else self.d)
    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])
    if warn and (self.bases==2).all():
        warnings.warn("It is more efficient to use DigitalNetB2 instead of DigitalNetAnyBases when all bases are 2")
    self.warn = warn

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 
53
54
55
56
57
58
59
60
61
62
63
64
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)

UML Specific