# Hyperuniformity

Functions designed to study the hyperuniformity of a stationary point process using estimation(s) of its structure factor StructureFactor.

Hyperuniformity diagnostics:

Note

Typical usage

1. Test the hyperuniformity using the statistical test multiscale_test() or the test of effective hyperuniformity effective_hyperuniformity().

2. If the hyperuniformity was approved, find the possible hyperuniformity class using hyperuniformity_class().

structure_factor.hyperuniformity.multiscale_test(point_pattern_list, estimator, subwindows_list, k_list, mean_poisson, m_list=None, proba_list=None, verbose=False, **kwargs)[source]

Compute the sample mean $$\bar{Z}$$ and the corresponding standard deviation $$\bar{\sigma}/\sqrt{N}$$ of the coupled sum estimator $$Z$$ of a point process using a list of $$N$$ PointPatterns and $$N$$ realizations from the random variable $$M$$.

The test rejects the hyperuniformity hypothesis if the confidence interval $$CI[\mathbb{E}[Z]]= [\bar{Z} - z \bar{\sigma}/\sqrt{N}, \bar{Z} + z \bar{\sigma}/\sqrt{N}]$$, doesn’t contain zero, and vice-versa. (See the multiscale hyperuniformity test in ).

Parameters
• point_pattern_list (list) – List of PointPattern (s). Each element of the list is an encapsulation of a realization of a point process, the observation window, and (optionally) the intensity of the point process. All PointPattern (s) should have the same window and intensity.

• estimator (str) – Choice of structure factor’s estimator. The parameters of the chosen estimator must be added as keyword arguments. The available estimators are “scattering_intensity”, “tapered_estimator”, “bartlett_isotropic_estimator”, and “quadrature_estimator_isotropic”. See StructureFactor.

• subwindows_list (list) – List of increasing cubic or ball-shaped AbstractSpatialWindow, typically, obtained using subwindows(). The shape of the windows depends on the choice of the estimator. Each element of point_pattern_list will be restricted to these windows to compute $$Z$$.

• k_list (list) – List of wavevectors (or wavenumbers) where the estimator is to be evaluated. Each element is associated with an element of subwindows_list. Typically, obtained using subwindows().

• mean_poisson (int) – Parameter of the Poisson r.v. $$M$$ used to compute $$Z$$. To use a different distribution of the r.v. $$M$$, set mean_poisson=None and specify m_list and proba_list corresponding to $$M$$.

• m_list (list, optional) – List of positive integer realizations of the r.v. $$M$$, used when mean_poisson=None. Each element of m_list is associated with an element of point_pattern_list to compute a realization of $$Z$$. Defaults to None.

• proba_list (list, optional) – List of $$\mathbb{P}(M \geq j)$$ used with m_list when mean_poisson=None. Should contains at least max(m_list) elements. Defaults to None.

• verbose (bool, optional) – If “True” and mean_poisson is not None, print the re-sampled values of $$M$$. Defaults to False.

Keyword Arguments

kwargs (dict) – Parameters of the chosen estimator of the structure factor. See StructureFactor.

Returns

• “mean_Z”: The sample mean of $$Z$$.

• ”std_mean_Z”: The sample standard deviation of $$Z$$ divided by the square root of the number of samples.

• ”Z”: The obtained values of $$Z$$.

• ”M”: The used values of $$M$$.

Return type

dict(float, float, list, list)

Example

import numpy as np
from structure_factor.point_processes import HomogeneousPoissonPointProcess
from structure_factor.spatial_windows import BoxWindow
from structure_factor.hyperuniformity import subwindows, multiscale_test

# PointPattern list
poisson = HomogeneousPoissonPointProcess(intensity=1 / np.pi)
N = 500
L = 160
window = BoxWindow([[-L / 2, L / 2]] * 2)
point_patterns = [
poisson.generate_point_pattern(window=window) for _ in range(N)
]
# subwindows and k
l_0 = 40
subwindows_list, k = subwindows(window, subwindows_type="BoxWindow", param_0=l_0)

# multiscale hyperuniformity test
mean_poisson = 90
summary = multiscale_test(point_patterns,
estimator="scattering_intensity", k_list=k,
subwindows_list=subwindows_list, mean_poisson=mean_poisson,
verbose=True)

import matplotlib.pyplot as plt
plt.hist(summary["Z"], color="grey", label="Z")
plt.axvline(summary["mean_Z"], color="b", linewidth=1, label=r"$\bar{Z}$")
plt.axvline(summary["mean_Z"] - 3*summary["std_mean_Z"], linestyle='dashed',
color="k", linewidth=1, label=r"$\bar{Z} \pm 3 \bar{\sigma}/ \sqrt{N}$")
plt.axvline(summary["mean_Z"] + 3*summary["std_mean_Z"], linestyle='dashed',
color="k", linewidth=1)
plt.xlabel(r"$Z$")
plt.title(r"Histogram of the $N$ obtained values of $Z$")
plt.legend()
plt.show()

Definition

Let $$\mathcal{X} \in \mathbb{R}^d$$ be a stationary point process of which we consider an increasing sequence of sets $$(\mathcal{X} \cap W_m)_{m \geq 1}$$, with $$(W_m)_m$$ centered box (or ball)-shaped windows s.t. $$W_s \subset W_r$$ for all $$0< s<r$$, and $$W_{\infty} = \mathbb{R}^d$$. We define the sequence of r.v. $$Y_m = 1\wedge \widehat{S}_m(\mathbf{k}_m^{\text{min}})$$, where $$\widehat{S}_m$$ is one of the positive, asymptotically unbiased estimators of the structure factor of $$\mathcal{X}$$ applied on the observation $$\mathcal{X} \cap W_m$$, and $$\mathbf{k}_m^{\text{min}}$$ is the minimum allowed wavevector associated with $$W_m$$.

Under some assumptions (, Section 4) $$\mathcal{X}$$ is hyperuniform iff $$\mathbb{E}[Z]=0$$. Where $$Z$$ is the coupled sum estimator of [RG15] defined by,

$Z = \sum_{j=1}^{M} \frac{Y_j - Y_{j-1}}{\mathbb{P}(M\geq j)},$

with $$M$$ an $$\mathbb{N}$$-valued random variable such that $$\mathbb{P}(M \geq j)>0$$ for all $$j$$, and $$Y_{0}=0$$.

Important

• If mean_poisson is not None, there is a step of accepting/rejecting while sampling from the r.v. $$M$$. If the biggest subwindow associated with $$M'$$ (obtained value of $$M$$) is larger than the father window, then $$M'$$ is rejected, and we resample from $$M$$. Typically, mean_poisson should be chosen s.t. the probability that the biggest subwindow is larger than the father window is small enough.

• The test is asymptotically valid, so it might fail in diagnosing hyperuniformity if the number of samples used to compute $$\bar{Z}$$ or $$\lambda$$ is too small.

Note

Typical usage

• The function subwindows() can be used to generate from a father window a list of subwindows and the associated allowed wavevectors/wavenumbers.

structure_factor.hyperuniformity.effective_hyperuniformity(k_norm, sf, k_norm_stop, std_sf=None, **kwargs)[source]

Evaluate the index $$H$$ of hyperuniformity of a point process using its structure factor sf. If $$H<10^{-3}$$ the corresponding point process is deemed effectively hyperuniform.

Parameters
• k_norm (numpy.ndarray) – Vector of wavenumbers (i.e. norms of wavevectors).

• sf (numpy.ndarray) – Evaluations of the structure factor, of the given point process, at k_norm.

• std_sf (numpy.ndarray, optional) – Standard deviations associated with sf. Defaults to None.

• k_norm_stop (float) – Threshold on k_norm. Used to find the numerator of $$H$$ by linear regression of sf up to the value associated with k_norm_stop.

Keyword Arguments

kwargs (dict) – Keyword arguments (except “sigma”) of scipy.scipy.optimize.curve_fit parameters.

Returns

• “H”: Value of the index $$H$$.

• ”s0_std”: Standard deviation of the numerator of $$H$$.

• ”fitted_line”: Line used to find the numerator of $$H$$.

• ”idx_first_peak”: Index of $$k_{peak}$$ if it exists.

Return type

dict(float, float, function, int)

Example

import matplotlib.pyplot as plt
import numpy as np
from structure_factor.hyperuniformity import effective_hyperuniformity
from structure_factor.point_processes import GinibrePointProcess
from structure_factor.structure_factor import StructureFactor
from structure_factor.tapered_estimators_isotropic import (
allowed_k_norm_bartlett_isotropic,
)

point_process = GinibrePointProcess()

sf = StructureFactor(point_pattern)
k_norm, sf_estimated = sf.bartlett_isotropic_estimator(k_norm)

summary = effective_hyperuniformity(k_norm, sf_estimated, k_norm_stop=0.2)
x = np.linspace(0, 1, 300)
sf_fitted_line = summary["fitted_line"](x)

sf_theoretical = point_process.structure_factor(k_norm)
fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(k_norm, sf_theoretical, "g", label=r"$S(\mathbf{k})$")
ax.plot(k_norm, sf_estimated, "b", marker=".", label="Approximated structure factor")
ax.plot(x, sf_fitted_line, "r--", label="Fitted line")

ax.annotate(
"H={}".format(summary["H"]),
xy=(0, 0),
xytext=(0.01, 0.1),
arrowprops=dict(facecolor="black", shrink=0.0001),
)
ax.legend()
ax.set_xlabel(r"$||\mathbf{k}||_2$")
ax.set_ylabel(r"$\mathsf{S}(\mathbf{k})$")

Definition

A stationary isotropic point process $$\mathcal{X} \subset \mathbb{R}^d$$, is said to be effectively hyperuniform if $$H \leq 10^{-3}$$ where $$H$$ is defined following [Tor18] (Section 11.1.6) and (supplementary Section 8) by,

$H = \frac{\hat{S}(\mathbf{0})}{S(\mathbf{k}_{peak})}\cdot$
• $$S$$ is the structure factor of $$\mathcal{X}$$,

• $$\hat{S}(\mathbf{0})$$ is a linear extrapolation of the structure factor at $$\mathbf{k}=\mathbf{0}$$,

• $$\mathbf{k}_{peak}$$ is the location of the first dominant peak value of $$S$$.

For more details, we refer to (Section 2.5).

Important

To compute $$\hat{S}(\mathbf{0})$$, a linear extrapolation with a least-square fit is used to fit a line on the values of sf associated with a subvector of k_norm. This subvector is obtained by truncating k_norm around k_norm_stop. For the choice of k_norm_stop, the trade-off is to remain close to zero while including enough data points to fit a line. In addition, std_sf will be considered while fitting the line if it’s not None.

structure_factor.hyperuniformity.hyperuniformity_class(k_norm, sf, k_norm_stop=1, std_sf=None, **kwargs)[source]

Fit a polynomial $$y = c \cdot x^{\alpha}$$ to sf around zero. $$\alpha$$ is used to specify the possible class of hyperuniformity of the associated point process (as described below).

Parameters
• k_norm (numpy.ndarray) – Vector of wavenumbers (i.e. norms of the wavevectors).

• sf (numpy.ndarray) – Evaluations of the structure factor, of the given point process, at k_norm.

• std (numpy.ndarray, optional) – Standard deviations associated to sf. Defaults to None.

• k_norm_stop (float, optional) – Threshold on k_norm. The subvector obtained from sf starting from zero up to the value associated with k_norm_stop is used to fit a polynomial and find $$\alpha$$. Defaults to 1.

Keyword Arguments

kwargs (dict) –

Keyword arguments (except “sigma”) of scipy.scipy.optimize.curve_fit parameters.

Returns

• “alpha”: Estimated value of $$\alpha$$.

• ”c”: Estimated value of $$c$$.

• ”fitted_poly”: Polynomial used to find $$\alpha$$.

Return type

dict(float, float, function)

Example

import matplotlib.pyplot as plt
from structure_factor.hyperuniformity import hyperuniformity_class
from structure_factor.point_processes import GinibrePointProcess
from structure_factor.structure_factor import StructureFactor
from structure_factor.tapered_estimators_isotropic import (
allowed_k_norm_bartlett_isotropic,
)

point_process = GinibrePointProcess()

sf = StructureFactor(point_pattern)
k_norm, sf_estimated = sf.bartlett_isotropic_estimator(k_norm)

summary = hyperuniformity_class(k_norm, sf_estimated, k_norm_stop=0.4)
sf_fitted_0 = summary["fitted_poly"](k_norm)

sf_theoretical = point_process.structure_factor(k_norm)

fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(k_norm, sf_theoretical, "g", label=r"$S(\mathbf{k})$")
ax.plot(k_norm, sf_estimated, "b", marker=".", label="Approximated structure factor")
ax.plot(k_norm, sf_fitted_0, "r--", label="Fitted line")

ax.annotate(
r"$\alpha$ ={}".format(summary["alpha"]),
xy=(0, 0),
xytext=(0.01, 0.1),
arrowprops=dict(facecolor="black", shrink=0.0001),
)
ax.legend()
ax.set_xlabel(r"$||\mathbf{k}||_2$")
ax.set_ylabel(r"$\mathsf{S}(\mathbf{k})$")


Definition

For a stationary hyperuniform point process $$\mathcal{X} \subset \mathbb{R}^d$$, if $$\vert S(\mathbf{k})\vert\sim c \Vert \mathbf{k} \Vert_2^\alpha$$ in the neighborhood of zero, then by [Cos21] (Section 4.1.) the value of $$\alpha$$ determines the hyperuniformity class of $$\mathcal{X}$$ as follows,

Class

$$\alpha$$

$$\mathbb{V}\text{ar}\left[\mathcal{X}(B(0, R)) \right]$$

I

$$> 1$$

$$\mathcal{O}(R^{d-1})$$

II

$$= 1$$

$$\mathcal{O}(R^{d-1}\log(R))$$

III

$$]0, 1[$$

$$\mathcal{O}(R^{d-\alpha})$$

For more details, we refer to , (Section 2.4).

Important

To compute $$\alpha$$, a polynomial is fitted on the values of sf associated with a subvector of k_norm. This subvector is obtained by truncating k_norm around k_norm_stop. For the choice of k_norm_stop, the trade-off is to remain close to zero while including enough data points to fit a polynomial. In addition, std_sf will be considered while fitting the polynomial if it’s not None.

structure_factor.hyperuniformity.bin_data(k_norm, sf, **params)[source]

Split k_norm into sub-intervals (or bins) and evaluate, over each sub-interval, the mean and the standard deviation of the corresponding values in sf.

Parameters
• k_norm (numpy.ndarray) – Vector of wavenumbers (i.e. norms of wavevectors).

• sf (numpy.ndarray) – Evaluations of the structure factor, of the given point process, at k_norm.

Keyword Arguments

params (dict) – Keyword arguments (except “x”, “values” and “statistic”) of scipy.stats.binned_statistic.

Returns

• k_norm: Centers of the bins.

• sf: Means of the structure factor over the bins.

• std_sf: Standard deviations of the structure factor over the bins.

Return type

tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray)

Example

import matplotlib.pyplot as plt
import numpy as np
import structure_factor.utils as utils
from structure_factor.hyperuniformity import bin_data
from structure_factor.point_processes import GinibrePointProcess
from structure_factor.spatial_windows import BoxWindow
from structure_factor.structure_factor import StructureFactor

point_process = GinibrePointProcess()

# Restrict point_pattern to window
window = BoxWindow([[-35, 35], [-35, 35]])
point_pattern_box = point_pattern.restrict_to_window(window)

# Estimate S
sf = StructureFactor(point_pattern_box)
x = np.linspace(0, 3, 80)
x = x[x != 0]
k = utils.meshgrid_to_column_matrix(np.meshgrid(x, x))
k, sf_estimated = sf.scattering_intensity(k)

# bin_data
k_norm = utils.norm(k)
k_norm_binned, sf_estimated_binned, _ = bin_data(k_norm, sf_estimated,bins=40)

fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(k_norm, sf_estimated, "b,", label="Before regularization", rasterized=True)
sf.plot_isotropic_estimator(
k_norm_binned,
sf_estimated_binned,
axis=ax,
color="m",
exact_sf=point_process.structure_factor,
label="After regularization",
)
ax.legend()


structure_factor.hyperuniformity.subwindows(window, subwindows_type='BoxWindow', param_0=None, param_max=None, params=None)[source]

Create a list of cubic (or ball)-shaped subwindows of a father window, with the associated minimum allowed wavevectors (or wavenumbers).

Parameters
• window (AbstractSpatialWindow) – Father window.

• subwindows_type (str, optional) – Type of the subwindows to be created. The available types are “BoxWindow” and “BallWindow”. The former for cubic and the latter for ball-shaped subwindows. Defaults to “BoxWindow”.

• param_0 (float, optional) – Parameter (lengthside/radius) of the first subwindow to be created. If not None, an increasing sequence of subwindows with parameters of unit increments is created. The biggest subwindow has parameter param_max if it’s not None, else, the maximum possible parameter. Defaults to None.

• param_max (float, optional) – Maximum subwindow parameter (lengthside/radius). Used when param_0 is not None. Defaults to None.

• params (list, optional) – List of parameters (lengthside/radius) of the output subwindows. For a list of parameters of unit increments, param_0 and param_max can be used instead. Defaults to None.

Returns

Return type

(list, list)

Example

from structure_factor.spatial_windows import BoxWindow
from structure_factor.hyperuniformity import subwindows

# Example 1:
L = 30
window = BoxWindow([[-L / 2, L / 2]] * 2)
# subwindows and k
l_0 = 10
subwindows_box, k = subwindows(window, subwindows_type="BoxWindow", param_0=l_0)

# Example 2:
subwindows_params=[1, 4, 7, 15]
subwindows_ball, k_norm = subwindows(window, subwindows_type="BallWindow", params=subwindows_params)

import matplotlib.pyplot as plt
fig, axis = plt.subplots(1, 2, figsize=(10,5))
axis.plot(0,0)
axis.plot(0,0)
for i, j in zip(subwindows_box, subwindows_ball):
i.plot(axis=axis)
j.plot(axis=axis)
plt.show()


Note

Typical usage