blue_sampler package

Submodules

blue_sampler.api module

Public API.

blue.sample_points is the main entry point of the package. It provides a high-level interface for generating large point sets on the periodic unit hypercube [0, 1)^D with sub-Poisson density fluctuations (so-called blue noise).

The package also exposes utilities to sample tessellations (2D only) and balanced clusters (arbitrary dimension). tessels or clusters are sampled with following balance property: uniform area repartition if no target is given or uniform atoms repartition if target atoms are given.

Clusters and tessellations can subsequently be converted into low-discrepancy point sets using `blue.from_geometry’ which internally solves a moment-matching problem

type blue_sampler.api.NDArray = ndarray[_AnyShape, dtype]

class blue_sampler.api.ProgressLogger(D, verbose)[source]

Bases: object

Hierarchical -based progress display for nested pipeline levels.

enter_level(N, D, N_ITER)[source]

Push a new recursion level and return its context.

Return type:: _LevelCtx

exit_level()[source]

Return type:: None

write(level, N, msg, newline=False)[source]

Return type:: None

blue_sampler.api.from_geometry(geometry, gtype, p=3, **kwargs)[source]

Convert tessels or clusters into low-discrepancy point set.

Each tessel or cluster is replaced by n points by solving a moment-matching problem (Levenberg-Marquardt) via momentum_fit.

Parameters:

geometry (NDArray) –
Batch of tessels or clusters to convert. - gtype=”polygons” : quadrilaterals of shape (N, 4, 2),

e.g. the quad output of sample_tessels.
- gtype=”clusters” : point sets of shape (N, K, D), e.g. the output of sample_clusters, or the atoms output of sample_tessels(…, return_atoms=True).
gtype (str) – Geometry type. “polygons” only supports D = 2.
p (int) – Maximum total moment order to match (centroid plus central moments up to order p).
**kwargs – Additional arguments passed to momentum_fit (e.g., n_restarts, random_state, tol).

Returns:

n points per tessel or cluster, matching its moments up to order p.

Return type:

NDArray

blue_sampler.api.im2points(image='anything.jpg', N=100000)[source]: simple wrapper for image stippling

blue_sampler.api.momentum_fit(distribution, distribution_type='clusters', p=3, weights=None, n_restarts=1, restart_tol=1e-30, n_iter=100, lambda0=0.01, tol=1e-30, random_state=None)[source]

Fit n points X_1, …, X_n whose centroid and central moments (q = 1, …, p-1) match those of distribution.

n is determined automatically (see required_n): it is not a parameter. When p == 2, only the centroid (q = 1) is requested; no central moments and no LM solve are needed, so the n = 1 target centroid is returned directly.

blue_sampler.api.plot(points, auto_zoom=False, max_points=30000, ax=None, return_fig=False, figsize=(8, 8), **scatter_kw)[source]

Scatter plot of a 2-D or 3-D point set.

If auto_zoom is set to True:: For large point sets the view is automatically zoomed so that at most max_points points are displayed.

Parameters:

points (NDArray) – Point coordinates, with D in {2, 3}. Higher-dimensional arrays are silently projected onto the first 3 axes.
auto_zoom (bool) – Whether to apply auto_zoom for large point sets.
max_points (int) – Maximum number of points to draw used for auto_zoom.
ax (Axes | None) – Existing axes to draw into. If None, a new figure is created.
return_fig (bool) – If True, returns (fig, ax). If False, displays the figure and returns None.
**scatter_kw – Extra keyword arguments forwarded to ax.scatter.

Return type:

tuple[Figure, Axes] | None

blue_sampler.api.sample_clusters(N=32768, D=2, targets=None, n_per_cluster=16)[source]

Recursively partition a point set into N balanced clusters.

At each recursion step, every cluster is split into two equal halves using a random median hyperplane. After log2(N) recursion levels, exactly N clusters are obtained.

Parameters:

N (int) – Number of final clusters. Must be a power of two.
D (int) – Ambient dimension.
targets (NDArray | None) – Initial atoms to clusterise. If omitted, a Sobol low-discrepancy sequence containing K = N * n_per_cluster atoms is generated automatically.
n_per_cluster (int) – Number of atoms per final cluster. Only used when targets is not provided.

Returns:

Collection of balanced clusters.

Return type:

NDArray

Notes

The recursive splitting procedure requires the total number of atoms K to be divisible by N.

blue_sampler.api.sample_points(N=32768, D=2, bruteforce=False, warmstart=None, n_iter=6, targets=None, verbose=1)[source]

Generate N stealthy (blue-noise) points in [0, 1)^D.

Parameters:

N (int) – Number of output points.
D (int) – Spatial dimension.
bruteforce (bool) – Use the bruteforce O(N^2) algorithm instead of the recursive one. Gives a better sample but is intractable for N >= 50_000. Automatically forced to True when N <= 2_000.
warmstart (str | NDArray | None) –
Initial point configuration. - None : default random/recursive initialisation. - “Sobol” : initialise with a Sobol low-discrepancy sequence

(requires scipy.stats.qmc.Sobol).
- ndarray : use given points as the starting configuration.
n_iter (int) – Number of solver iterations. Each iteration runs 10 gradient steps plus a structural gridification step (neighbor lookup). More iterations = better quality but slower.
targets (NDArray | None) – Atoms describing a target density. can also be a path to an image, e.g. targets = “zebra.jpg”
verbose (int) – 0 = silent, 1 = live progress.

Returns:

points – The sampled point coordinates in [0, 1)^D.

Return type:

NDArray

blue_sampler.api.sample_tessels(N=32768, D=2, targets=None, return_atoms=False)[source]

Recursively split the unit square into N random quadrilaterals.

If targets is None, splits are computed to achieve equal areas. If targets is provided, splits are computed to achieve a median separation of the atoms.

Parameters:

N (int) – Number of final tessels. Must be a power of 2.
targets (NDArray | None) – Coordinates of atoms to split, in the [0, 1)^2 unit box. K must be a multiple of N. The more targets provided, the better the approximation, but the slower the computation (K/N should be at least 100 for a decent tessellation). A typical use case is adaptive tessellation, e.g. with targets being i.i.d. points sampled from a target density.

Return type:

NDArray | tuple[NDArray, NDArray]

Returns:

quad (ndarray of shape (N, 4, 2)) – A tessellation composed of N quadrilaterals with equal area or equal number of atoms.
atoms (ndarray of shape (N, K/N, 2)) – The input target atoms, redistributed among their final quadrilateral. Only returned if targets was provided and return_atoms is set to True.

Notes

Only supports 2D geometry and a power-of-two number of tessels (N).

blue_sampler.api.tile(x, repeat, flatoutput=True)[source]

Tile points on the unit torus to cover [0, 1)^D periodically.

Each of the repeat**D copies of x is rescaled by 1/repeat and shifted to its own sub-cube, so that the copies together pave the unit torus again. For example in 2D with repeat=2: tile (0, 0) holds x/2, tile (1, 1) holds x/2 + 0.5, etc.

Parameters:

x (NDArray) – Points in [0, 1)^D (unit torus).
repeat (int) – Number of repetitions per axis. The output therefore contains Nfinal = N * repeat**D points.
flatoutput (bool) – If True, reshape the output to (Nfinal, D). If False, keep the tile structure as leading axes.

Returns:

Tiled version of x, periodized over [0, 1)^D.

Return type:

NDArray

blue_sampler.datasets module

blue_sampler.datasets.generate_dataset(size=(1000, 2), dataset='barbara', plot_data=True, list_method=True)[source]

blue_sampler.datasets.list_methods()[source]: Return available sampling methods.

blue_sampler.datasets.samplepoints(method='gaussian', size=(1000000, 2), plot_points=False)[source]

Generate point clouds using various sampling strategies.

Parameters:

method (str) –
Sampling method. Supported values:

General methods (any dimension D): - “square” : uniform in [-1, 1]^D - “ball” : uniform in the unit ball - “ring” : Gaussian with radial offset - “gaussian” : standard normal distribution - “spiky” : L^alpha ball (alpha < 1) - “holy” : hypercube with spherical holes

Dataset-based methods (fixed dimensions): - “barbara”, “lena”, “france”, “germany” : 2D datasets - “everest” : 3D dataset (elevation-based sampling)
size (tuple of int) – (N, D) where N is the number of points and D the dimension. Note: D must match the dataset dimension for dataset-based methods.
plot_points (bool, optional) – If True, display the generated points.

Returns:

Array of shape (N, D) containing sampled points.

Return type:

numpy.ndarray

blue_sampler.math module

Low-level mathematical helpers

blue_sampler.math.clean_grad(x)[source]

Replace NaN gradient contributions (fictive points) with 0.

Return type:: Array

blue_sampler.math.drop_symmetric(directions)[source]

Keep only one representative from each direction pair {v, -v}.

The canonical representative is the one whose first non-zero component is positive.

Parameters:: directions (ndarray)
Return type:: ndarray

blue_sampler.math.grid_shape(N, D)[source]

Smallest D-hypercube grid that contains at least N points.

Return type:

tuple[tuple[int, ...], int, tuple[int, ...]]

Returns:

IJK (shape tuple e.g. (32, 32) for D=2)
total (total number of grid slots (I^D))
axes (tuple(range(D)))

blue_sampler.math.integers_in_half_ball(radius, D)[source]

Return all non-zero integer lattice vectors inside a sphere of radius, keeping only one vector per direction pair.

Parameters:

radius (float)
D (int)

Return type:

ndarray

blue_sampler.math.prepare_points(x, N_asked, IJK, D)[source]

Pad N_asked real points to fill the I^D grid.

Fictive slots receive a NaN coordinate so gradients ignore them.

Parameters:

x (ndarray | None)
N_asked (int)
IJK (tuple[int, ...])
D (int)

Return type:

Array

blue_sampler.math.prepare_wave_vectors(Ks)[source]

Build JAX arrays for the spectral gradient.

Parameters:

Ks (ndarray)

Return type:

tuple[Array, Array]

Returns:

K_w (complex array of shape (M, D) — phase multipliers)
K_ (complex array of shape (M, D) — normalised duals)

blue_sampler.math.random_rotations(x, batch_size, Dout, Din)[source]

blue_sampler.math.sample_wave_vectors(kmed, kmax, D, n_high)[source]

Return type:: ndarray

blue_sampler.math.simplex(D)[source]

Vertices of a regular simplex centred at the origin in R^D.

Return type:: ndarray

blue_sampler.math.structure_factor(points, resolution=2000)[source]

Estimate the radial structure factor S(k) via scattering intensity.

Parameters:

points (ndarray)
resolution (int)

Return type:

tuple[ndarray, ndarray]

Returns:

k ((M,) float array — exact wave-vector magnitudes.)
S ((M,) float array — exact S(k) values.)

Note

beyond a radius fixed to capture 1/4 of the “resolution” budget, ALL allowed wavevectors are sampled to get maximal precision on low frequencies. Remaining budget is shared evenly accross all pertinent frequency scales.

blue_sampler.math.torus_delta(delta)[source]

Shortest signed displacement on the unit torus.

Return type:: Array

blue_sampler.math.torus_wrap(x)[source]

Wrap coordinates into [0, 1)^D.

Return type:: Array

blue_sampler.progress module

Just some progress loging staff used when verbose is set to 1 in the solver

class blue_sampler.progress.ProgressLogger(D, verbose)[source]

Bases: object

Hierarchical -based progress display for nested pipeline levels.

enter_level(N, D, N_ITER)[source]

Push a new recursion level and return its context.

Return type:: _LevelCtx

exit_level()[source]

Return type:: None

write(level, N, msg, newline=False)[source]

Return type:: None

blue_sampler.viz module

Visualisation helpers: plot display a 2-D or 3-D point set plot_structure_factor display the structure factor estimated through

scattering intensity

plot_tessels display a tessellation at up to 12 recursion depths plot_clusters display a cluster partition at up to 12 recursion depths

type blue_sampler.viz.NDArray = ndarray[_AnyShape, dtype]

blue_sampler.viz.back_merge_clusters(clusters)[source]

Inverse operation of a recursion level.

` Merge sibling clusters back together following the ordering convention used by ``fair_random_split.

blue_sampler.viz.back_merge_tessels(quads, eps=1e-07)[source]: Invert one step of the tesselation tree by pairing brother quads and identifying the split orientation using cyclic vertex permutation. 100% vectorized.

blue_sampler.viz.plot(points, auto_zoom=False, max_points=30000, ax=None, return_fig=False, figsize=(8, 8), **scatter_kw)[source]

Scatter plot of a 2-D or 3-D point set.

If auto_zoom is set to True:: For large point sets the view is automatically zoomed so that at most max_points points are displayed.

Parameters:

points (NDArray) – Point coordinates, with D in {2, 3}. Higher-dimensional arrays are silently projected onto the first 3 axes.
auto_zoom (bool) – Whether to apply auto_zoom for large point sets.
max_points (int) – Maximum number of points to draw used for auto_zoom.
ax (Axes | None) – Existing axes to draw into. If None, a new figure is created.
return_fig (bool) – If True, returns (fig, ax). If False, displays the figure and returns None.
**scatter_kw – Extra keyword arguments forwarded to ax.scatter.

Return type:

tuple[Figure, Axes] | None

blue_sampler.viz.plot_clusters(clusters, return_fig=False)[source]

Display a cluster partition across up to 12 recursion steps.

Return type:: tuple[Figure, NDArray] | None

blue_sampler.viz.plot_structure_factor(points, resolution=2000, smoothed=True, ax=None, return_fig=False, **plot_kw)[source]

Log-log plot of the radial structure factor S(k).

S(k) is estimated using standard scattering intensity.

Parameters:

points (NDArray) – Point coordinates in [0, 1)^D.
resolution (int) – Number of sampled wave vectors used to estimate sf.
smoothed (bool) – If True, apply a local log-log average. If False, plot raw values.
ax (Axes | None) – Existing axes to draw into. If None, a new figure is created.
return_fig (bool) – If True, returns (fig, ax).
**plot_kw – Extra keyword arguments forwarded to ax.loglog.

Return type:

tuple[Figure, Axes] | None

blue_sampler.viz.plot_tessels(tessels, return_fig=False)[source]

Display a tessellation across up to 12 recursion steps.

Return type:: tuple[Figure, NDArray] | None

blue_sampler.viz.show_clusters(ax, clusters)[source]

Return type:: Axes

blue_sampler.viz.show_polygons(ax, tessels)[source]

Return type:: Axes

blue_sampler.warm_start module

Module contents

blue_sampler

Generate stealthy point patterns — low-discrepancy, spectrally isotropic samples on the unit torus [0, 1)^D.

Quick start

>>> import blue_sampler as blue
>>> x = blue.sample(N=10_000, D=2)          # (10000, 2) array
>>> blue.plot(x)
>>> blue.plot_structure_factor(x)

blue_sampler.from_geometry(geometry, gtype, p=3, **kwargs)[source]

Convert tessels or clusters into low-discrepancy point set.

Each tessel or cluster is replaced by n points by solving a moment-matching problem (Levenberg-Marquardt) via momentum_fit.

Parameters:

geometry (NDArray) –
Batch of tessels or clusters to convert. - gtype=”polygons” : quadrilaterals of shape (N, 4, 2),

e.g. the quad output of sample_tessels.
- gtype=”clusters” : point sets of shape (N, K, D), e.g. the output of sample_clusters, or the atoms output of sample_tessels(…, return_atoms=True).
gtype (str) – Geometry type. “polygons” only supports D = 2.
p (int) – Maximum total moment order to match (centroid plus central moments up to order p).
**kwargs – Additional arguments passed to momentum_fit (e.g., n_restarts, random_state, tol).

Returns:

n points per tessel or cluster, matching its moments up to order p.

Return type:

NDArray

blue_sampler.generate_dataset(size=(1000, 2), dataset='barbara', plot_data=True, list_method=True)[source]

blue_sampler.im2points(image='anything.jpg', N=100000)[source]: simple wrapper for image stippling

blue_sampler.plot(points, auto_zoom=False, max_points=30000, ax=None, return_fig=False, figsize=(8, 8), **scatter_kw)[source]

Scatter plot of a 2-D or 3-D point set.

If auto_zoom is set to True:: For large point sets the view is automatically zoomed so that at most max_points points are displayed.

Parameters:

points (NDArray) – Point coordinates, with D in {2, 3}. Higher-dimensional arrays are silently projected onto the first 3 axes.
auto_zoom (bool) – Whether to apply auto_zoom for large point sets.
max_points (int) – Maximum number of points to draw used for auto_zoom.
ax (Axes | None) – Existing axes to draw into. If None, a new figure is created.
return_fig (bool) – If True, returns (fig, ax). If False, displays the figure and returns None.
**scatter_kw – Extra keyword arguments forwarded to ax.scatter.

Return type:

tuple[Figure, Axes] | None

blue_sampler.plot_clusters(clusters, return_fig=False)[source]

Display a cluster partition across up to 12 recursion steps.

Return type:: tuple[Figure, NDArray] | None

blue_sampler.plot_structure_factor(points, resolution=2000, smoothed=True, ax=None, return_fig=False, **plot_kw)[source]

Log-log plot of the radial structure factor S(k).

S(k) is estimated using standard scattering intensity.

Parameters:

points (NDArray) – Point coordinates in [0, 1)^D.
resolution (int) – Number of sampled wave vectors used to estimate sf.
smoothed (bool) – If True, apply a local log-log average. If False, plot raw values.
ax (Axes | None) – Existing axes to draw into. If None, a new figure is created.
return_fig (bool) – If True, returns (fig, ax).
**plot_kw – Extra keyword arguments forwarded to ax.loglog.

Return type:

tuple[Figure, Axes] | None

blue_sampler.plot_tessels(tessels, return_fig=False)[source]

Display a tessellation across up to 12 recursion steps.

Return type:: tuple[Figure, NDArray] | None

blue_sampler.sample_clusters(N=32768, D=2, targets=None, n_per_cluster=16)[source]

Recursively partition a point set into N balanced clusters.

At each recursion step, every cluster is split into two equal halves using a random median hyperplane. After log2(N) recursion levels, exactly N clusters are obtained.

Parameters:

N (int) – Number of final clusters. Must be a power of two.
D (int) – Ambient dimension.
targets (Optional[NDArray]) – Initial atoms to clusterise. If omitted, a Sobol low-discrepancy sequence containing K = N * n_per_cluster atoms is generated automatically.
n_per_cluster (int) – Number of atoms per final cluster. Only used when targets is not provided.

Returns:

Collection of balanced clusters.

Return type:

NDArray

Notes

The recursive splitting procedure requires the total number of atoms K to be divisible by N.

blue_sampler.sample_points(N=32768, D=2, bruteforce=False, warmstart=None, n_iter=6, targets=None, verbose=1)[source]

Generate N stealthy (blue-noise) points in [0, 1)^D.

Parameters:

N (int) – Number of output points.
D (int) – Spatial dimension.
bruteforce (bool) – Use the bruteforce O(N^2) algorithm instead of the recursive one. Gives a better sample but is intractable for N >= 50_000. Automatically forced to True when N <= 2_000.
warmstart (Union[str, NDArray, None]) –
Initial point configuration. - None : default random/recursive initialisation. - “Sobol” : initialise with a Sobol low-discrepancy sequence

(requires scipy.stats.qmc.Sobol).
- ndarray : use given points as the starting configuration.
n_iter (int) – Number of solver iterations. Each iteration runs 10 gradient steps plus a structural gridification step (neighbor lookup). More iterations = better quality but slower.
targets (Optional[NDArray]) – Atoms describing a target density. can also be a path to an image, e.g. targets = “zebra.jpg”
verbose (int) – 0 = silent, 1 = live progress.

Returns:

points – The sampled point coordinates in [0, 1)^D.

Return type:

NDArray

blue_sampler.sample_tessels(N=32768, D=2, targets=None, return_atoms=False)[source]

Recursively split the unit square into N random quadrilaterals.

If targets is None, splits are computed to achieve equal areas. If targets is provided, splits are computed to achieve a median separation of the atoms.

Parameters:

N (int) – Number of final tessels. Must be a power of 2.
targets (Optional[NDArray]) – Coordinates of atoms to split, in the [0, 1)^2 unit box. K must be a multiple of N. The more targets provided, the better the approximation, but the slower the computation (K/N should be at least 100 for a decent tessellation). A typical use case is adaptive tessellation, e.g. with targets being i.i.d. points sampled from a target density.

Return type:

Union[NDArray, tuple[NDArray, NDArray]]

Returns:

quad (ndarray of shape (N, 4, 2)) – A tessellation composed of N quadrilaterals with equal area or equal number of atoms.
atoms (ndarray of shape (N, K/N, 2)) – The input target atoms, redistributed among their final quadrilateral. Only returned if targets was provided and return_atoms is set to True.

Notes

Only supports 2D geometry and a power-of-two number of tessels (N).

blue_sampler.structure_factor(points, resolution=2000)[source]

Estimate the radial structure factor S(k) via scattering intensity.

Parameters:

points (ndarray)
resolution (int)

Return type:

tuple[ndarray, ndarray]

Returns:

k ((M,) float array — exact wave-vector magnitudes.)
S ((M,) float array — exact S(k) values.)

Note

beyond a radius fixed to capture 1/4 of the “resolution” budget, ALL allowed wavevectors are sampled to get maximal precision on low frequencies. Remaining budget is shared evenly accross all pertinent frequency scales.

blue_sampler.tile(x, repeat, flatoutput=True)[source]

Tile points on the unit torus to cover [0, 1)^D periodically.

Each of the repeat**D copies of x is rescaled by 1/repeat and shifted to its own sub-cube, so that the copies together pave the unit torus again. For example in 2D with repeat=2: tile (0, 0) holds x/2, tile (1, 1) holds x/2 + 0.5, etc.

Parameters:

x (NDArray) – Points in [0, 1)^D (unit torus).
repeat (int) – Number of repetitions per axis. The output therefore contains Nfinal = N * repeat**D points.
flatoutput (bool) – If True, reshape the output to (Nfinal, D). If False, keep the tile structure as leading axes.

Returns:

Tiled version of x, periodized over [0, 1)^D.

Return type:

NDArray