Sobol & Halton Generators

pointpats/random.py

@@ -25,6 +25,7 @@
 
 import numpy
 from scipy.spatial import cKDTree
+from scipy.stats import qmc
 
 from .geometry import area as _area
 from .geometry import bbox as _bbox
@@ -700,3 +701,288 @@ def strauss(
             )
 
     return result.squeeze()
+
+
+def sobol(hull, intensity=None, size=None, scramble=True, seed=None):
+    """
+    Simulate a point pattern using a Sobol low-discrepancy sequence.
+
+    Parameters
+    ----------
+    hull : A geometry-like object
+        This encodes the "space" in which to simulate the pattern.
+        All points will lie within this hull. Supported values are:
+
+        - a bounding box encoded as ``numpy.array([xmin, ymin, xmax, ymax])``
+        - an (N,2) array of points for which the bounding box will be computed & used
+        - a shapely polygon/multipolygon
+        - a scipy convex hull
+
+    intensity : float
+        The number of observations per unit area in the hull to use. If
+        provided, then ``size`` must be an integer describing the number of
+        replications to generate.
+
+    size : tuple or int
+        A tuple of ``(n_observations, n_replications)``, where the first
+        number is the number of points to simulate in each replication and
+        the second number is the number of replications. If an integer is
+        provided and ``intensity`` is None, one replication is generated.
+        If an integer is provided and ``intensity`` is also supplied, the
+        integer specifies the number of replications and the number of
+        observations is computed from the intensity.
+
+    scramble : bool, default True
+        If True, apply Owen scrambling to the Sobol sequence. Scrambling
+        preserves the low-discrepancy properties while allowing independent
+        realizations to be generated using different seeds.
+
+    seed : int or None, optional
+        Seed used to initialize the scrambled Sobol sequence. Ignored when
+        ``scramble=False``.
+
+    Returns
+    -------
+    numpy.ndarray
+        Either an ``(n_replications, n_observations, 2)`` array or an
+        ``(n_observations, 2)`` array containing the simulated realizations.
+
+    Examples
+    --------
+    Generate 100 Sobol points in the unit square.
+
+    >>> points = sobol(numpy.array([0, 0, 1, 1]), size=100, seed=123)
+    >>> points.shape
+    (100, 2)
+
+    Generate 3 independent replications.
+
+    >>> points = sobol(numpy.array([0, 0, 1, 1]), size=(100, 3), seed=123)
+    >>> points.shape
+    (3, 100, 2)
+    """
+    if isinstance(hull, numpy.ndarray):
+        hull = hull if hull.shape == (4,) else _prepare_hull(hull)
+
+    n_observations, n_simulations, intensity = parse_size_and_intensity(
+        hull, intensity=intensity, size=size
+    )
+
+    result = numpy.empty((n_simulations, n_observations, 2))
+
+    bbox = _bbox(hull)
+
+    for i_replication in range(n_simulations):
+        sbl = qmc.Sobol(
+            d=2,
+            scramble=scramble,
+            seed=None if seed is None else seed + i_replication,
+        )
+
+        accepted = []
+
+        while len(accepted) < n_observations:
+            remaining = n_observations - len(accepted)
+
+            m = int(numpy.ceil(numpy.log2(max(remaining, 1))))
+            candidates = sbl.random_base2(m=m)
+
+            xs = bbox[0] + candidates[:, 0] * (bbox[2] - bbox[0])
+            ys = bbox[1] + candidates[:, 1] * (bbox[3] - bbox[1])
+
+            for x, y in zip(xs, ys, strict=True):
+                if _contains(hull, x, y):
+                    accepted.append((x, y))
+                if len(accepted) == n_observations:
+                    break
+
+        result[i_replication] = numpy.asarray(accepted)
+
+    return result.squeeze()
+
+
+def halton(hull, intensity=None, size=None, scramble=True, seed=None):
+    """
+    Simulate a point pattern using a Halton low-discrepancy sequence.
+
+    Parameters
+    ----------
+    hull : A geometry-like object
+        This encodes the "space" in which to simulate the pattern.
+        All points will lie within this hull. Supported values are:
+
+        - a bounding box encoded as ``numpy.array([xmin, ymin, xmax, ymax])``
+        - an (N,2) array of points for which the bounding box will be computed & used
+        - a shapely polygon/multipolygon
+        - a scipy convex hull
+
+    intensity : float
+        The number of observations per unit area in the hull to use. If
+        provided, then ``size`` must be an integer describing the number of
+        replications to generate.
+
+    size : tuple or int
+        A tuple of ``(n_observations, n_replications)``, where the first
+        number is the number of points to simulate in each replication and
+        the second number is the number of replications. If an integer is
+        provided and ``intensity`` is None, one replication is generated.
+        If an integer is provided and ``intensity`` is also supplied, the
+        integer specifies the number of replications and the number of
+        observations is computed from the intensity.
+
+    scramble : bool, default True
+        If True, apply scrambling to the Halton sequence. Scrambling reduces
+        correlation artifacts and allows independent realizations to be
+        generated using different seeds.
+
+    seed : int or None, optional
+        Seed used to initialize the scrambled Halton sequence. Ignored when
+        ``scramble=False``.
+
+    Returns
+    -------
+    numpy.ndarray
+        Either an ``(n_replications, n_observations, 2)`` array or an
+        ``(n_observations, 2)`` array containing the simulated realizations.
+
+    Examples
+    --------
+    Generate 100 Halton points in the unit square.
+
+    >>> points = halton(numpy.array([0, 0, 1, 1]), size=100, seed=123)
+    >>> points.shape
+    (100, 2)
+
+    Generate 3 independent replications.
+
+    >>> points = halton(numpy.array([0, 0, 1, 1]), size=(100, 3), seed=123)
+    >>> points.shape
+    (3, 100, 2)
+    """
+    if isinstance(hull, numpy.ndarray):
+        hull = hull if hull.shape == (4,) else _prepare_hull(hull)
+
+    n_observations, n_simulations, intensity = parse_size_and_intensity(
+        hull, intensity=intensity, size=size
+    )
+
+    result = numpy.empty((n_simulations, n_observations, 2))
+
+    bbox = _bbox(hull)
+
+    for i_replication in range(n_simulations):
+        hltn = qmc.Halton(
+            d=2,
+            scramble=scramble,
+            seed=None if seed is None else seed + i_replication,
+        )
+
+        accepted = []
+
+        while len(accepted) < n_observations:
+            remaining = n_observations - len(accepted)
+
+            candidates = hltn.random(remaining)
+
+            xs = bbox[0] + candidates[:, 0] * (bbox[2] - bbox[0])
+            ys = bbox[1] + candidates[:, 1] * (bbox[3] - bbox[1])
+
+            for x, y in zip(xs, ys, strict=True):
+                if _contains(hull, x, y):
+                    accepted.append((x, y))
+
+        result[i_replication] = numpy.asarray(accepted)
+
+    return result.squeeze()
+
+
+def r2(hull, intensity=None, size=None, seed=None):
+    """
+    Simulate a point pattern using the Roberts R2 low-discrepancy sequence.
+
+    Parameters
+    ----------
+    hull : A geometry-like object
+        This encodes the "space" in which to simulate the pattern.
+        All points will lie within this hull. Supported values are:
+
+        - a bounding box encoded as ``numpy.array([xmin, ymin, xmax, ymax])``
+        - an (N,2) array of points for which the bounding box will be computed & used
+        - a shapely polygon/multipolygon
+        - a scipy convex hull
+
+    intensity : float
+        The number of observations per unit area in the hull to use. If
+        provided, then ``size`` must be an integer describing the number of
+        replications to generate.
+
+    size : tuple or int
+        A tuple of ``(n_observations, n_replications)``, where the first
+        number is the number of points to simulate in each replication and
+        the second number is the number of replications. If an integer is
+        provided and ``intensity`` is None, one replication is generated.
+        If an integer is provided and ``intensity`` is also supplied, the
+        integer specifies the number of replications and the number of
+        observations is computed from the intensity.
+
+    seed : int or None, optional
+        Seed used to initialize the R2 sequence. If None, a random seed
+        is used.
+
+    Returns
+    -------
+    numpy.ndarray
+        Either an ``(n_replications, n_observations, 2)`` array or an
+        ``(n_observations, 2)`` array containing the simulated realizations.
+
+    Examples
+    --------
+    Generate 100 R2 points in the unit square.
+
+    >>> points = r2(numpy.array([0, 0, 1, 1]), size=100, seed=123)
+    >>> points.shape
+    (100, 2)
+
+    Generate 3 independent replications.
+
+    >>> points = r2(numpy.array([0, 0, 1, 1]), size=(100, 3), seed=123)
+    >>> points.shape
+    (3, 100, 2)
+    """
+    if isinstance(hull, numpy.ndarray):
+        hull = hull if hull.shape == (4,) else _prepare_hull(hull)
+
+    n_observations, n_simulations, intensity = parse_size_and_intensity(
+        hull, intensity=intensity, size=size
+    )
+
+    result = numpy.empty((n_simulations, n_observations, 2))
+
+    bbox = _bbox(hull)
+
+    phi2 = 1.324717957244746
+    alpha = numpy.array([1.0 / phi2, 1.0 / phi2**2])
+
+    for i_replication in range(n_simulations):
+        rng_seed = None if seed is None else seed + i_replication
+        offset = float(numpy.random.default_rng(rng_seed).random())
+        accepted = []
+        n_generated = 0
+
+        while len(accepted) < n_observations:
+            remaining = n_observations - len(accepted)
+
+            indices = numpy.arange(n_generated + 1, n_generated + remaining + 1)
+            candidates = (offset + numpy.outer(indices, alpha)) % 1.0
+            n_generated += remaining
+
+            xs = bbox[0] + candidates[:, 0] * (bbox[2] - bbox[0])
+            ys = bbox[1] + candidates[:, 1] * (bbox[3] - bbox[1])
+
+            for x, y in zip(xs, ys, strict=True):
+                if _contains(hull, x, y):
+                    accepted.append((x, y))
+
+        result[i_replication] = numpy.asarray(accepted)
+
+    return result.squeeze()

pointpats/tests/test_random.py

@@ -0,0 +1,42 @@
+import numpy as np
+
+from pointpats.random import halton, sobol
+
+
+# -----------------------------------------------------------------------------
+# Sobol tests
+# -----------------------------------------------------------------------------
+def test_sobol_shape():
+    hull = np.array([0, 0, 1, 1])
+    points = sobol(hull, size=10)
+    assert points.shape == (10, 2)
+
+def test_sobol_multiple_replications():
+    hull = np.array([0, 0, 1, 1])
+    points = sobol(hull, size=(10, 4))
+    assert points.shape == (4, 10, 2)
+
+def test_sobol_reproducible():
+    hull = np.array([0, 0, 1, 1])
+    p1 = sobol(hull, size=100, seed=123)
+    p2 = sobol(hull, size=100, seed=123)
+    np.testing.assert_allclose(p1, p2)
+
+# -----------------------------------------------------------------------------
+# Halton tests
+# -----------------------------------------------------------------------------
+def test_halton_shape():
+    hull = np.array([0, 0, 1, 1])
+    points = halton(hull, size=10)
+    assert points.shape == (10, 2)
+
+def test_halton_multiple_replications():
+    hull = np.array([0, 0, 1, 1])
+    points = halton(hull, size=(10, 4))
+    assert points.shape == (4, 10, 2)
+
+def test_halton_reproducible():
+    hull = np.array([0, 0, 1, 1])
+    p1 = halton(hull, size=100, seed=123)
+    p2 = halton(hull, size=100, seed=123)
+    np.testing.assert_allclose(p1, p2)