From e573a85ceb91596b79d8d1f817f5e08d5a8dd995 Mon Sep 17 00:00:00 2001 From: R Virinchi Date: Thu, 18 Jun 2026 10:15:12 +0530 Subject: [PATCH 1/5] Added sobol and halton qmc generators --- pointpats/random.py | 96 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) diff --git a/pointpats/random.py b/pointpats/random.py index 7ff01051..2963f47a 100644 --- a/pointpats/random.py +++ b/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,98 @@ def strauss( ) return result.squeeze() + +def sobol( + hull, + intensity=None, + size=None, + scramble=True, + seed=None): + """ + todo + """ + 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 + + 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): + 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): + """ + todo + """ + 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 + + 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): + if _contains(hull, x, y): + accepted.append((x, y)) + + result[i_replication] = numpy.asarray(accepted) + + return result.squeeze() From 68c1196d0a5fddebcffeb59791498151d6d9bace Mon Sep 17 00:00:00 2001 From: R Virinchi Date: Thu, 18 Jun 2026 10:15:54 +0530 Subject: [PATCH 2/5] new file for testing random.py ruff fixes added strict parameter to zip fixed formatting fixed --- pointpats/random.py | 132 ++++++++++++++++++++++++++++++--- pointpats/tests/test_random.py | 42 +++++++++++ 2 files changed, 162 insertions(+), 12 deletions(-) create mode 100644 pointpats/tests/test_random.py diff --git a/pointpats/random.py b/pointpats/random.py index 2963f47a..e8dd10bb 100644 --- a/pointpats/random.py +++ b/pointpats/random.py @@ -709,18 +709,72 @@ def sobol( scramble=True, seed=None): """ - todo + 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([0, 0, 1, 1], size=100, seed=123) + >>> points.shape + (100, 2) + + Generate 3 independent replications. + + >>> points = sobol([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 + bbox = _bbox(hull) for i_replication in range(n_simulations): sbl = qmc.Sobol( @@ -732,7 +786,7 @@ def sobol( accepted = [] while len(accepted) < n_observations: - + remaining = n_observations - len(accepted) m = int(numpy.ceil(numpy.log2(max(remaining, 1)))) @@ -741,14 +795,14 @@ def sobol( 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): + 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( @@ -758,18 +812,72 @@ def halton( scramble=True, seed=None): """ - todo + 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([0, 0, 1, 1], size=100, seed=123) + >>> points.shape + (100, 2) + + Generate 3 independent replications. + + >>> points = halton([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 + bbox = _bbox(hull) for i_replication in range(n_simulations): hltn = qmc.Halton( @@ -777,7 +885,7 @@ def halton( scramble=scramble, seed=None if seed is None else seed + i_replication, ) - + accepted = [] while len(accepted) < n_observations: @@ -789,7 +897,7 @@ def halton( 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): + for x, y in zip(xs, ys, strict=True): if _contains(hull, x, y): accepted.append((x, y)) diff --git a/pointpats/tests/test_random.py b/pointpats/tests/test_random.py new file mode 100644 index 00000000..5e8a4411 --- /dev/null +++ b/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) From 4fdefe1883f8469bb90af4ea22e6c28d227ac37a Mon Sep 17 00:00:00 2001 From: R Virinchi Date: Thu, 18 Jun 2026 23:55:38 +0530 Subject: [PATCH 3/5] fixed indentation --- pointpats/random.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/pointpats/random.py b/pointpats/random.py index e8dd10bb..00c29a9c 100644 --- a/pointpats/random.py +++ b/pointpats/random.py @@ -886,20 +886,20 @@ def halton( seed=None if seed is None else seed + i_replication, ) - accepted = [] + accepted = [] - while len(accepted) < n_observations: + while len(accepted) < n_observations: - remaining = n_observations - len(accepted) + remaining = n_observations - len(accepted) - candidates = hltn.random(remaining) + candidates = hltn.random(remaining) - xs = bbox[0] + candidates[:, 0] * (bbox[2] - bbox[0]) - ys = bbox[1] + candidates[:, 1] * (bbox[3] - bbox[1]) + 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)) + for x, y in zip(xs, ys, strict=True): + if _contains(hull, x, y): + accepted.append((x, y)) result[i_replication] = numpy.asarray(accepted) From 3a60d960f91bdbc2f2220eb726d152598b7f7c92 Mon Sep 17 00:00:00 2001 From: R Virinchi Date: Wed, 1 Jul 2026 00:07:14 +0530 Subject: [PATCH 4/5] added implementation of r2 generator --- pointpats/random.py | 102 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 98 insertions(+), 4 deletions(-) diff --git a/pointpats/random.py b/pointpats/random.py index 00c29a9c..d1a5bdf7 100644 --- a/pointpats/random.py +++ b/pointpats/random.py @@ -702,6 +702,7 @@ def strauss( return result.squeeze() + def sobol( hull, intensity=None, @@ -755,13 +756,13 @@ def sobol( -------- Generate 100 Sobol points in the unit square. - >>> points = sobol([0, 0, 1, 1], size=100, seed=123) + >>> points = sobol(numpy.array([0, 0, 1, 1]), size=100, seed=123) >>> points.shape (100, 2) Generate 3 independent replications. - >>> points = sobol([0, 0, 1, 1], size=(100, 3), seed=123) + >>> points = sobol(numpy.array([0, 0, 1, 1]), size=(100, 3), seed=123) >>> points.shape (3, 100, 2) """ @@ -805,6 +806,7 @@ def sobol( return result.squeeze() + def halton( hull, intensity=None, @@ -858,13 +860,13 @@ def halton( -------- Generate 100 Halton points in the unit square. - >>> points = halton([0, 0, 1, 1], size=100, seed=123) + >>> points = halton(numpy.array([0, 0, 1, 1]), size=100, seed=123) >>> points.shape (100, 2) Generate 3 independent replications. - >>> points = halton([0, 0, 1, 1], size=(100, 3), seed=123) + >>> points = halton(numpy.array([0, 0, 1, 1]), size=(100, 3), seed=123) >>> points.shape (3, 100, 2) """ @@ -904,3 +906,95 @@ def halton( 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() From 5a3a7b244e3f19556e52ed9be0139d5a37063bc8 Mon Sep 17 00:00:00 2001 From: R Virinchi Date: Wed, 1 Jul 2026 00:10:14 +0530 Subject: [PATCH 5/5] fixed ruff errors --- pointpats/random.py | 18 +++--------------- 1 file changed, 3 insertions(+), 15 deletions(-) diff --git a/pointpats/random.py b/pointpats/random.py index d1a5bdf7..fb532c7f 100644 --- a/pointpats/random.py +++ b/pointpats/random.py @@ -703,12 +703,7 @@ def strauss( return result.squeeze() -def sobol( - hull, - intensity=None, - size=None, - scramble=True, - seed=None): +def sobol(hull, intensity=None, size=None, scramble=True, seed=None): """ Simulate a point pattern using a Sobol low-discrepancy sequence. @@ -787,7 +782,6 @@ def sobol( accepted = [] while len(accepted) < n_observations: - remaining = n_observations - len(accepted) m = int(numpy.ceil(numpy.log2(max(remaining, 1)))) @@ -807,12 +801,7 @@ def sobol( return result.squeeze() -def halton( - hull, - intensity=None, - size=None, - scramble=True, - seed=None): +def halton(hull, intensity=None, size=None, scramble=True, seed=None): """ Simulate a point pattern using a Halton low-discrepancy sequence. @@ -891,7 +880,6 @@ def halton( accepted = [] while len(accepted) < n_observations: - remaining = n_observations - len(accepted) candidates = hltn.random(remaining) @@ -983,7 +971,7 @@ def r2(hull, intensity=None, size=None, seed=None): 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