📥 click here to download as a Jupyter notebook

0. Installation

Don't forget to clone the github repository:

git clone https://github.com/For-a-few-DPPs-more/hyperuniform-samplers.git 
cd hyperuniform-samplers

Then install numpy and h5py:

%pip install numpy h5py

I. Load the datasets from Zenodo

The database has two families, selected by the ball keyword:

ball	Spectral exclusion region	Use for
l2	Euclidean ball (isotropic)	blue noise, rendering, spatial sampling
wd	wrap-around discrepancy set (axis-oriented)	QMC integration, low effective dimension

Both use χ = 0.43 and guarantee S(k) <= 1e-10 on their exclusion region (numerical zero in float32).

Open a new terminal (outside this notebook) and run from the root of the project:

Option 1: one (or several) specific dimensions and ball (--ball accepts l2, wd or all)

python src/load_data.py --dim 2 3 4 --ball l2
python src/load_data.py --dim 2 3 4 --ball wd

Option 2: the full database
⚠️ about 800 MB, may take several minutes.

python src/load_data.py --dim all --ball all

Files are saved in data/<ball>/D.../, you can follow the progress from there.
Once the download is finished, come back to this notebook and continue.

II. Create a new point cloud

A new randomized point cloud is generated from the downloaded one by exploiting the symmetries of the distribution:

• random shifts
• random flips
• random axis swaps

When needed (see VI. Take Home), projection and tiling are applied first, taking a point cloud from shape (N, D) to (N × 2^i, D − j).

import os
import numpy as np
from src.new_data import new_data

os.makedirs("new_data", exist_ok=True)
X = new_data(N=1024, D=2, ball="l2", seed=42)
np.save("new_data/N1024_D2_l2_seed42.npy", X)

To run new_data from the terminal instead:

python src/new_data.py --dim 2 --npoints 1024 --ball l2 --seed 42
# other family: --ball wd

III. Analysis

This requires the blue-sampler python package (%pip install blue-sampler).

import blue_sampler as blue
from src.new_data import new_data

for ball in ("l2", "wd"):
    X = new_data(N=2048, D=2, ball=ball, seed=42)
    blue.plot(X)
    blue.plot_structure_factor(X)

IV. Spectral landscape (2D only)

import numpy as np
import matplotlib.pyplot as plt
from src.new_data import new_data


def plot2D_structure_factor(pts, wv, ax=None, cmap="gray", figsize=(7, 6)):
    """
    Plot the 2D structure factor S(k) at selected wavevectors.

    Parameters
    ----------
    pts : array-like, shape (N, 2)
        Positions of the points in real space.
    wv : array-like, shape (M, 2)
        Wavevectors at which to evaluate S(k).
    ax : matplotlib.axes.Axes, optional
        Existing axes.
    cmap : str
        Colormap used for S(k).
    figsize : tuple
        Figure size if ax is None.

    Returns
    -------
    fig, ax, sf : figure, axes, ndarray
        Figure, axes and computed structure factor.
    """
    pts = np.asarray(pts, dtype=float)
    wv = np.asarray(wv, dtype=float)
    wv = wv[~(wv == 0).all(axis=-1)]
    # S(k) = |sum_j exp(2i pi k.r_j)|^2 / N
    phase = 2 * np.pi * wv @ pts.T
    sf = np.abs(np.exp(1j * phase).sum(axis=1)) ** 2 / len(pts)

    if ax is None:
        fig, ax = plt.subplots(figsize=figsize)
    else:
        fig = ax.figure

    wv = np.concatenate((wv, -wv), axis=0)
    sf = np.concatenate((sf, sf), axis=0)
    sf = np.log(sf).clip(min=-10, max=3)

    sc = ax.scatter(
        wv[:, 0], wv[:, 1],
        c=sf,
        cmap=cmap,
        s=300 * (1000 / len(wv)),
        linewidths=0,
        edgecolors="none",
    )

    fig.colorbar(sc, ax=ax)
    ax.set_title("2D Structure Factor (log scale)")
    ax.set_xticks([])
    ax.set_yticks([])

    return fig, ax, sf


N = 1024
D = 2

Nsq = int(np.sqrt(N))
wv_x = np.arange(Nsq)
wv_y = np.arange(-Nsq + 1, Nsq)
wv = np.stack(np.meshgrid(wv_x, wv_y, indexing="ij"), axis=-1).reshape(-1, 2)

for ball in ("l2", "wd"):
    pts = new_data(N, D, ball=ball, seed=42)
    fig, ax, _ = plot2D_structure_factor(pts, wv)
    ax.set_title(f"2D Structure Factor (log scale), ball = {ball}")

V. Benchmark with other qmc methods

This requires pip install scipy
From the terminal, run: python src/benchmark_data.py
This will compute scipy.stats.qmc.discrepancy() for several samplers.

VI. Take Home

Download the datasets (from the root of the repository), either manually from Zenodo (DOI 10.5281/zenodo.23244509)
or by running, in the terminal:

# everything: D = 2 to 16, plus 23 and 32, both families (l2, wd)
python src/load_data.py --dim all --ball all
# specific dimensions, specific ball (faster)
python src/load_data.py --dim 2 3 4 --ball l2
# the wd family, e.g. for QMC integration
python src/load_data.py --dim 2 3 4 --ball wd

Build a new randomised (N, D) point cloud from the downloaded ones:

• from Python:

from src.new_data import new_data
X = new_data(N=1024, D=2, ball="l2", seed=42, data_dir="data")  # (N, D) np.float32 array
X_wd = new_data(N=1024, D=2, ball="wd", seed=42)                # wd family

• or from the terminal:

python src/new_data.py --dim 2 --npoints 1024 --ball l2 --seed 42
python src/new_data.py --dim 2 --npoints 1024 --ball wd --seed 42

The ball parameter (l2 or wd) selects the dataset. Both families use χ = 0.43 and guarantee S(k) <= 1e-10 on their exclusion region, which is the numerical zero in float32:

• l2: Euclidean-ball exclusion region (isotropic, good blue-noise properties). Default choice.
• wd: exclusion region defined by the wrap-around discrepancy (axis-oriented), recommended for QMC integration of integrands with low effective dimension.

new_data.py performs, in this order:

1. Projection: if D is not stored (17 <= D <= 22 or 24 <= D <= 31), a dataset with D' > D is used and projected.
2. Tiling: if N >= 2^18 (about 260k), the largest stored cloud (N = 2^17) is tiled, see TILING.md.
3. Randomisation (with the given seed): random shift x_new = (x - s) mod 1, random flips x_new[:, d] = 1 - x_new[:, d],
   and random axis swaps.

Supported: any N power of 2 with 2^10 <= N <= 2^24, any 2 <= D <= 32, ball = {l2, wd}.

There is no hard limit on tiling, but we recommend not going far beyond N = 2^24:
stealthiness degrades proportionally to N_t / N (the S(k) <= 1e-10 guarantee holds for the native, untiled samples).