Source code for particles.qmc

# -*- coding: utf-8 -*-

"""QMC and RQMC sequences.

This module is a simple wrapper for LowDiscrepancy.f, a fortran piece of code
that implements the Sobol' and Halton sequence. (The same fortran code
is used the `randtoolbox
<>` package
in ``R``.)

To use this module, you must first compile LowDiscrepancy.f, as explained in
the installation notes.

import numpy as np
import warnings

qmc_warning = """
Module lowdiscrepancy does not seem to be installed (see INSTALL
notes). You will not be able to run SQMC (quasi-Monte Carlo version
of SMC).

sobol_warning =  """
lowdiscrepancy.sobol(%i, %i, %i, %i, 1, 0) generated points outside (0, 1)
    from particles import lowdiscrepancy
except ImportError:

max_int_32 = np.iinfo(np.int32).max

[docs]def sobol(N, dim, scrambled=1): """ Sobol sequence. Parameters ---------- N : int length of sequence dim: int dimension scrambled: int which scrambling method to use: + 0: no scrambling + 1: Owen's scrambling + 2: Faure-Tezuka + 3: Owen + Faure-Tezuka Returns ------- (N, dim) numpy array. Notes ----- One of the argument of the underlying Fortran function is the seed for the random generator used for scrambling. We simply generate this seed uniformly (between 0 and 2^32 - 1). There is a very small number of seeds that generate points that are == 0 (This has been reported to the maintainer of randtoolbox). When this happens, we generate a warning and start over (i.e. we re-generate another random seed, and compute a new scrambled Sobol point set. """ while True: seed = np.random.randint(max_int_32) out = lowdiscrepancy.sobol(N, dim, scrambled, seed, 1, 0) if np.all(out < 1.) and np.all(out > 0.): return out else: warnings.warn(sobol_warning % (N, dim, scrambled, seed))
[docs]def halton(N, dim): """ Halton sequence. Component i of the sequence consists of a Van der Corput sequence in base b_i, where b_i is the i-th prime number. Parameters ---------- N : int length of sequence dim: int dimension Returns ------- (N, dim) numpy array. """ return lowdiscrepancy.halton(N, dim, 1, 0)