"""Collection of point processes at related properties, e.g., intensity, pair correlation function, structure factor.
- :py:class:`~structure_factor.point_processes.HomogeneousPoissonPointProcess`: The homogeneous Poisson point process.
- :py:class:`~structure_factor.point_processes.ThomasPointProcess`: The Thomas point process.
- :py:class:`~structure_factor.point_processes.GinibrePointProcess`: The Ginibre point process.
- :py:func:`~structure_factor.point_processes.mutual_nearest_neighbor_matching`: The matching process of :cite:`KlaLasYog20`.
"""
import numpy as np
import scipy.linalg as la
from scipy.spatial import KDTree
from structure_factor.point_pattern import PointPattern
from structure_factor.spatial_windows import (
AbstractSpatialWindow,
BallWindow,
BoxWindow,
)
from structure_factor.utils import get_random_number_generator
[docs]class HomogeneousPoissonPointProcess(object):
"""`Homogeneous Poisson point process <https://en.wikipedia.org/wiki/Poisson_point_process#Spatial_Poisson_point_process>`_.
.. todo::
list attributes
"""
[docs] def __init__(self, intensity=1.0):
"""Create a `homogeneous Poisson point process <https://en.wikipedia.org/wiki/Poisson_point_process#Spatial_Poisson_point_process>`_ with prescribed (positive) ``intensity`` parameter.
Args:
intensity (float, optional): intensity of the homogeneous Poisson point process. Defaults to 1.0.
"""
if not intensity > 0:
raise TypeError("intensity argument must be 2positive")
self._intensity = float(intensity)
@property
def intensity(self):
r"""Return the intensity :math:`\rho_1(r) = \rho` of the homogeneous Poisson point process.
Returns:
float: Constant intensity.
"""
return self._intensity
[docs] @staticmethod
def pair_correlation_function(r=None):
r"""Evaluate the pair correlation function :math:`g(r)=1` of the homogeneous Poisson process.
Args:
r (numpy.ndarray): Array of size :math:`n \times d`, where :math:`d` is the ambient dimension and :math:`n` the number of points where the pair correlation function is evaluated. Since the homogeneous Poisson process is isotropic, a vector of size :math:`n` corresponding to the norm of the :math:`n` points can also be passed. Defaults to None.
Returns:
float or numpy.ndarray: ``1.0`` if ``r=None``, otherwise a vector of size :math:`n` with entries equal to ``1.0``.
"""
val = 1.0
if r is None:
return val
assert r.ndim <= 2
return np.full((r.shape[0], 1), val)
[docs] @staticmethod
def structure_factor(k=None):
r"""Evaluate the structure factor :math:`S(k)=1` of the homogeneous Poisson process.
Args:
k (numpy.ndarray): Array of size :math:`n \times d`, where :math:`d` is the ambient dimension and :math:`n` the number of points where the structure factor is evaluated. Since the homogeneous Poisson process is isotropic, a vector of size :math:`n` corresponding to the norm of the :math:`n` points can also be passed. Defaults to None.
Returns:
float or numpy.ndarray: ``1.0`` if ``k=None``, otherwise a vector of size :math:`n` with entries equal to ``1.0``.
"""
val = 1.0
if k is None:
return val
assert k.ndim <= 2
return np.full(k.shape[0], val)
[docs] def generate_sample(self, window, seed=None):
r"""Generate an exact sample of the point process restricted to the :math:`d` dimensional `window`.
Args:
window (:py:class:`~structure_factor.spatial_windows.AbstractSpatialWindow`): :math:`d`-dimensional observation window where the points will be generated.
Returns:
numpy.ndarray: Array of size :math:`n \times d`, where :math:`n` is the number of points forming the sample.
Example:
.. plot:: code/point_processes/poisson_sample.py
:include-source: True
:caption:
:alt: code/point_processes/poisson_sample.py
:align: center
.. seealso::
- :py:mod:`~structure_factor.spatial_windows`
"""
if not isinstance(window, AbstractSpatialWindow):
raise TypeError("window argument must be an AbstractSpatialWindow")
rng = get_random_number_generator(seed)
rho = self.intensity
nb_points = rng.poisson(rho * window.volume)
return window.rand(nb_points, seed=rng)
[docs] def generate_point_pattern(self, window, seed=None):
"""Generate a :py:class:`~structure_factor.point_pattern.PointPattern` of the point process.
Args:
window (:py:class:`~structure_factor.spatial_windows.AbstractSpatialWindow`): Observation window.
seed (int, optional): Seed to initialize the points generator. Defaults to None.
Returns:
:py:class:`~structure_factor.point_pattern.PointPattern`:Object of type :py:class:`~structure_factor.point_pattern.PointPattern` containing a realization of the point process, the observation window, and (optionally) the intensity of the point process (see :py:class:`~structure_factor.point_pattern.PointPattern`).
Example:
.. plot:: code/point_processes/poisson_pp.py
:include-source: True
:caption:
:alt: code/point_processes/poisson_pp.py
:align: center
.. seealso::
- :py:mod:`~structure_factor.spatial_windows`
- :py:class:`~structure_factor.point_pattern.PointPattern`
"""
points = self.generate_sample(window=window, seed=seed)
point_pattern = PointPattern(
points=points, window=window, intensity=self.intensity
)
return point_pattern
[docs]class ThomasPointProcess:
"""Homogeneous Thomas point process with Gaussian clusters.
.. todo::
list attributes
"""
[docs] def __init__(self, kappa, mu, sigma):
"""Create a homogeneous Thomas point process.
Args:
kappa (float): Mean number of clusters per unit volume. Intensity of the parent Poisson point process.
mu (float): Mean number of points per clusters. Mean of the child Gaussian distribution.
sigma (float): Standard deviation of the child Gaussian distribution.
"""
self.kappa = kappa
self.mu = mu
self.sigma = sigma
self._intensity = float(kappa * mu)
@property
def intensity(self):
r"""Return the intensity :math:`\rho_1(r) = \kappa \mu` of the Thomas point process.
Returns:
float: Constant intensity.
"""
return self._intensity
[docs] def pair_correlation_function(self, r_norm, d=2):
r"""Evaluate the pair correlation function of the Thomas point process.
.. math::
g(r)
= 1
+ \kappa (4 \pi \sigma^2)^{d / 2} \exp(-\frac{1}{4\sigma^2} \|r\|^2)
Args:
r_norm (numpy.ndarray): Vector of size :math:`n` corresponding to the norm of the :math:`n` points where the pair correlation function is evaluated.
d (int, optional): Ambient dimension. Defaults to 2.
Returns:
numpy.ndarray: Vector of size :math:`n` containing the evaluation of the pair correlation function.
"""
k = self.kappa
s2 = self.sigma ** 2
pcf = np.exp(r_norm ** 2 / (-4 * s2))
pcf /= k * (4 * np.pi * s2) ** (d / 2)
pcf += 1.0
return pcf
[docs] def structure_factor(self, k_norm):
r"""Evaluate the structure factor of the Thomas point process.
.. math::
S(k) = 1 + \mu \exp(- \sigma^2 * \|k\|^2).
Args:
k_norm (numpy.ndarray): Vector of size :math:`n` corresponding to the norm of the :math:`n` points where the structure factor is evaluated.
Returns:
numpy.ndarray: Vector of size :math:`n` containing the evaluation of the structure factor.
"""
mu = self.mu
s2 = self.sigma ** 2
return 1.0 + mu * np.exp(-s2 * k_norm ** 2)
[docs] def generate_sample(self, window, seed=None):
r"""Generate an exact sample from the corresponding :py:class:`~structure_factor.thomas_process.ThomasPointProcess` restricted to the :math:`d` dimensional `window`.
Args:
window (:py:class:`~structure_factor.spatial_windows.AbstractSpatialWindow`): :math:`d`-dimensional observation window where the points will be generated.
Returns:
numpy.ndarray: Array of size :math:`n \times d`, where :math:`n` is the number of points forming the sample.
Example:
.. plot:: code/point_processes/thomas_sample.py
:include-source: True
:caption:
:alt: code/point_processes/thomas_sample.py
:align: center
.. seealso::
- :py:mod:`~structure_factor.spatial_windows`
"""
if not isinstance(window, AbstractSpatialWindow):
raise TypeError("window argument must be an AbstractSpatialWindow")
rng = get_random_number_generator(seed)
tol = 6 * self.sigma
if isinstance(window, BoxWindow):
extended_bounds = window.bounds + np.array([-tol, +tol])
extended_window = BoxWindow(extended_bounds)
elif isinstance(window, BallWindow):
exented_radius = window.radius + tol
extended_window = BallWindow(window.center, exented_radius)
pp = HomogeneousPoissonPointProcess(self.kappa)
centers = pp.generate_sample(extended_window, seed=rng)
n_per_cluster = np.random.poisson(self.mu, size=len(centers))
d = window.dimension
s = self.sigma
points = np.vstack(
[rng.normal(c, s, (n, d)) for (c, n) in zip(centers, n_per_cluster)]
)
mask = window.indicator_function(points)
return points[mask]
[docs] def generate_point_pattern(self, window, seed=None):
"""Generate a :py:class:`~structure_factor.point_pattern.PointPattern` of the point process.
Args:
window (:py:class:`~structure_factor.spatial_windows.AbstractSpatialWindow`): Observation window.
seed (int, optional): Seed to initialize the points generator. Defaults to None.
Returns:
:py:class:`~structure_factor.point_pattern.PointPattern`:Object of type :py:class:`~structure_factor.point_pattern.PointPattern` containing a realization of the point process, the observation window, and (optionally) the intensity of the point process (see :py:class:`~structure_factor.point_pattern.PointPattern`).
Example:
.. plot:: code/point_processes/thomas_pp.py
:include-source: True
:caption:
:alt: code/point_processes/thomas_pp.py
:align: center
.. seealso::
- :py:mod:`~structure_factor.spatial_windows`
- :py:class:`~structure_factor.point_pattern.PointPattern`
"""
points = self.generate_sample(window=window, seed=seed)
point_pattern = PointPattern(
points=points, window=window, intensity=self.intensity
)
return point_pattern
[docs]class GinibrePointProcess(object):
"""Ginibre point process corresponds to the complex eigenvalues of a standard complex Gaussian matrix.
.. todo::
list attributes
"""
[docs] def __init__(self):
self._intensity = 1.0 / np.pi
@property
def intensity(self):
r"""Return the intensity :math:`\rho_1(r) = \frac{1}{\pi}` of the Ginibre point process.
Returns:
float: Constant intensity.
"""
return self._intensity
[docs] @staticmethod
def pair_correlation_function(r_norm):
r"""Evaluate the pair correlation function of the Ginibre point process.
.. math::
g(r) = 1 - \exp(- \|r\|^2)
Args:
r_norm (numpy.ndarray): Vector of size :math:`n` corresponding to the norm of the :math:`n` points where the pair correlation function is evaluated.
Returns:
numpy.ndarray: Vector of size :math:`n` containing the evaluation of the pair correlation function.
"""
return 1.0 - np.exp(-(r_norm ** 2))
[docs] @staticmethod
def structure_factor(k_norm):
r"""Evaluate the structure factor of the Ginibre point process.
.. math::
S(k) = 1 - \exp(- \frac{1}{4} \|k\|^2).
Args:
k_norm (numpy.ndarray): Vector of size :math:`n` corresponding to the norm of the :math:`n` points where the structure factor is evaluated.
Returns:
numpy.ndarray: Vector of size :math:`n` containing the evaluation of the structure factor.
"""
return 1.0 - np.exp(-0.25 * (k_norm ** 2))
[docs] def generate_sample(self, window, n=None, seed=None):
r"""Generate a sample of the Ginibre point process made of ``n`` points, inside the observation ``window``.
This is done by computing the eigenvalues of a random matrix :math:`G` of size :math:`n \times n`, filled with i.i.d. standard complex Gaussian variables, i.e.,
.. math::
G_{ij} = \frac{1}{\sqrt{2}} (X_{ij} + \mathbf{i} Y_{ij}).
Args:
window (:py:class:`~structure_factor.spatial_windows.BallWindow`): :math:`2`-dimensional centered ball window where the points will be generated.
n (int, optional): Number of points of the output sample. If ``n`` is None (default), it is set to the integer part of :math:`\rho |W| = \frac{1}{\pi} |W|`. Defaults to None.
Returns:
numpy.ndarray: Array of size :math:`n \times 2`, representing the :math:`n` points forming the sample.
Example:
.. plot:: code/point_processes/ginibre_sample.py
:include-source: True
:caption:
:alt: code/point_processes/ginibre_sample.py
:align: center
.. seealso::
- :py:class:`~structure_factor.spatial_windows.BallWindow`
"""
if not isinstance(window, BallWindow):
raise ValueError("The window should be a 2-d centered BallWindow.")
if window.dimension != 2:
raise ValueError("The window should be a 2-d window.")
if not np.all(np.equal(window.center, 0.0)):
raise ValueError("The window should be a centered window.")
if n is None:
n = int(window.volume * self.intensity)
assert isinstance(n, int)
rng = get_random_number_generator(seed)
A = np.zeros((n, n), dtype=complex)
A.real = rng.standard_normal((n, n))
A.imag = rng.standard_normal((n, n))
eigvals = la.eigvals(A) / np.sqrt(2.0)
points = np.vstack((eigvals.real, eigvals.imag)).T
mask = window.indicator_function(points)
return points[mask]
[docs] def generate_point_pattern(self, window, n=None, seed=None):
r"""Generate a :math:`2`-dimensional :py:class:`~structure_factor.point_pattern.PointPattern` of the point process, with a centered :py:class:`~structure_factor.spatial_windows.BallWindow`.
Args:
window (:py:class:`~structure_factor.spatial_windows.AbstractSpatialWindow`): :math:`2`-dimensional observation centered :py:class:`~structure_factor.spatial_windows.BallWindow`.
n (int, optional): Number of points of the output sample. If ``n`` is None (default), it is set to the integer part of :math:`\rho |W| = \frac{1}{\pi} |W|`. Defaults to None.
seed (int, optional): Seed to initialize the points generator. Defaults to None.
Returns:
:py:class:`~structure_factor.point_pattern.PointPattern`:Object of type :py:class:`~structure_factor.point_pattern.PointPattern` containing a realization of the point process, the observation window, and (optionally) the intensity of the point process (see :py:class:`~structure_factor.point_pattern.PointPattern`).
Example:
.. plot:: code/point_processes/ginibre_pp.py
:include-source: True
:caption:
:alt: code/point_processes/ginibre_pp.py
:align: center
.. seealso::
- :py:class:`~structure_factor.spatial_windows.BallWindow`
- :py:class:`~structure_factor.point_pattern.PointPattern`
"""
points = self.generate_sample(window=window, n=n, seed=seed)
point_pattern = PointPattern(
points=points, window=window, intensity=self.intensity
)
return point_pattern
[docs]def mutual_nearest_neighbor_matching(X, Y, **KDTree_params):
r"""Match the set of points ``X`` with a subset of points from ``Y`` based on mutual nearest neighbor matching :cite:`KlaLasYog20`. It is assumed that :math:`|X| \leq |Y|` and that each point in ``X``, resp. ``Y``, can have only one nearest neighbor in ``Y``, resp. ``X``.
The matching routine involves successive 1-nearest neighbor sweeps performed by `scipy.spatial.KDTree <https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.KDTree.html>`_ with the euclidean distance.
Args:
X (numpy.ndarray): Array of size (m, d) collecting points to be matched with a subset of points from ``Y``.
Y (numpy.ndarray): Array of size (n, d) of points satisfying :math:`m \leq n`.
Keyword Args:
see (documentation): keyword arguments of `scipy.spatial.KDTree <https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.KDTree.html>`_.
.. note::
The ``boxsize`` keyword argument can be used **only** when points belong to a box :math:`\prod_{i=1}^{d} [0, L_i)` (upper boundary excluded). It allows to consider periodic boundaries, i.e., the toroidal distance is used for searching for nearest neighbors.
Returns:
numpy.ndarray: vector of indices ``matches`` such that ``X[i]`` is matched to ``Y[matches[i]]``.
Example:
.. plot:: code/point_processes/kly_matching.py
:include-source: True
:caption: KLY matching
:alt: code/point_processes/kly_matching.py
:align: center
"""
if not (X.ndim == Y.ndim == 2):
raise ValueError(
"X and Y must be 2d numpy arrays with respective size (m, d) and (n, d), where d is the ambient dimension."
)
if X.shape[0] > Y.shape[0]:
raise ValueError(
"The sets of points represented by X and Y must satisfy |X| <= |Y|."
)
m, n = X.shape[0], Y.shape[0]
idx_X_unmatched = np.arange(m, dtype=int)
idx_Y_unmatched = np.arange(n, dtype=int)
matches = np.zeros(m, dtype=int)
for _ in range(m): # at most |X| nearest neighbor sweeps are performed
X_ = X[idx_X_unmatched]
Y_ = Y[idx_Y_unmatched]
knn = KDTree(Y_, **KDTree_params)
X_to_Y = knn.query(X_, k=1, p=2)[1] # p=2, i.e., euclidean distance
knn = KDTree(X_, **KDTree_params)
Y_to_X = knn.query(Y_, k=1, p=2)[1]
identity = range(len(idx_X_unmatched))
mask_X = np.equal(Y_to_X[X_to_Y], identity)
matches[idx_X_unmatched[mask_X]] = idx_Y_unmatched[X_to_Y[mask_X]]
if np.all(mask_X): # all points from X got matched
break
idx_X_unmatched = idx_X_unmatched[~mask_X]
mask_Y = np.full(len(idx_Y_unmatched), True, dtype=bool)
mask_Y[X_to_Y[mask_X]] = False
idx_Y_unmatched = idx_Y_unmatched[mask_Y]
return matches