Source code for structure_factor.tapered_estimators_isotropic

"""Collection of functions used to implement a fast version of Bartlett's isotropic estimator."""


import numpy as np
import scipy.special as sc
from numba import njit
from scipy.spatial.distance import pdist

from structure_factor.spatial_windows import BallWindow, UnitBallWindow


#! care about case k=0
[docs]def bartlett_estimator(k_norm, point_pattern): r"""Compute an estimation of the structure factor of a stationary isotropic point process from one realization encapsulated in ``point_pattern``, evaluated at ``k_norm``. Args: point_pattern (:py:class:`~structure_factor.point_pattern.PointPattern`): Realization of the underlying point process. k_norm (numpy.ndarray, optional): Array of size :math:`n` corresponding to the wavenumbers where the estimator is to be evaluated. If None (default) and the observation window ``point_pattern.window`` is a :py:class:`~structure_factor.spatial_windows.BallWindow` and the ambient dimension is even, the estimator will be evaluated on the corresponding set of allowed wavenumbers returned by :py:func:`~structure_factor.tapered_estimators_isotropic.allowed_k_norm`. Defaults to None. n_allowed_k_norm (int, optional): Number of allowed wavenumbers to be used when ``k_norm`` is None. See :py:func:`~structure_factor.tapered_estimators_isotropic.allowed_k_norm`. Defaults to 60. Returns: tuple(numpy.ndarray, numpy.ndarray): - k_norm: Wavenumber(s) at which the estimator has been evaluated. - estimation: Evaluation(s) of the estimator at ``k_norm``. Example: .. literalinclude:: code/isotropic_estimator/bartlett_estimator.py :language: python .. proof:definition:: The Bartlett's isotropic estimator :math:`\widehat{S}_{\mathrm{BI}}` is constructed from a realization :math:`\{\mathbf{x}_i\}_{i=1}^N` of the point process :math:`\mathcal{X}` observed in a ball window :math:`W=B(\mathbf{0}, R)`. .. math:: \widehat{S}_{\mathrm{BI}}(k) = 1 + \frac{ (2\pi)^{d/2} }{\rho |W| \omega_{d-1}} \sum_{ \substack{j, q =1 \\ j\neq q } }^{N} \frac{1}{(k \|\mathbf{x}_j - \mathbf{x}_q\|_2)^{d/2 - 1}} J_{d/2 - 1}(k \|\mathbf{x}_j - \mathbf{x}_q\|_2). For more details, we refer to :cite:`HGBLR:22`, (Section 3.2). .. seealso:: - :py:class:`~structure_factor.spatial_windows.BallWindow` - :py:meth:`~structure_factor.point_pattern.PointPattern.restrict_to_window` - :py:meth:`~structure_factor.structure_factor.StructureFactor.tapered_estimator_isotropic` - :py:func:`~structure_factor.tapered_estimators_isotropic.allowed_k_norm` """ if not isinstance(point_pattern.window, BallWindow): raise TypeError( "The window must be an instance of BallWindow. Hint: use point_pattern.restrict_to_window." ) window = point_pattern.window d = window.dimension X = np.atleast_2d(point_pattern.points) norm_xi_xj = pdist(X, metric="euclidean") k_norm = k_norm.astype(float) sf_estimated = np.zeros_like(k_norm) order = float(d / 2 - 1) bartlett_estimator_njit(sf_estimated, k_norm, norm_xi_xj, order) surface = UnitBallWindow(np.zeros(d)).surface volume = window.volume rho = point_pattern.intensity sf_estimated *= (2.0 * np.pi) ** (d / 2) / (surface * volume * rho) sf_estimated += 1.0 return sf_estimated
[docs]def allowed_k_norm_bartlett_isotropic(dimension, radius, nb_values=60): r"""Allowed wavenumbers of the Bartlett isotropic estimator, for a ``d``-dimensional point process observed in a ball window with radius ``radius``. Only for even values of ``d``. .. warning:: This method is only available when the ambient dimension is even. Args: dimension (int): Dimension (even number) of the underlying space. radius (float): Radius of the observation window. nb_values (int): Number of required allowed wavenumbers. Defaults to 60. Returns: numpy.ndarray: Vector of size ``nb_values`` containing the allowed wavenumbers. Example: .. literalinclude:: code/isotropic_estimator/allowed_k_norm.py :language: python .. testoutput:: [0.1915853 0.35077933 0.50867341 0.6661846 0.8235315 0.98079293] .. proof:definition:: The allowed wavenumbers of a realization from a point process :math:`\mathcal{X}` of :math:`\mathbb{R}^d` observed in a ball window :math:`W=B(\mathbf{0}, R)` correspond to .. math:: \left\{ \frac{x}{R} \in \mathbb{R} \text{ s.t. } J_{d/2}(x)=0 \right\} .. seealso:: - :py:func:`~structure_factor.tapered_estimators_isotropic.bartlett_estimator` """ d_2, mod_ = divmod(dimension, 2) is_even_dimension = mod_ == 0 if not is_even_dimension: # dimension is not even raise ValueError( "Allowed wavenumber could be used only when the dimension of the space `d` is an even number (i.e., d/2 is an integer)." ) return sc.jn_zeros(d_2, nb_values) / radius
@njit("double(double, double)") def bessel1_njit(order, x): if order == 0.0: return sc.j0(x) if order == 1.0: return sc.j1(x) return sc.jv(order, x) # todo use cuda.jit to decrease computational time https://towardsdatascience.com/python-performance-and-gpus-1be860ffd58d#:~:text=x%5B%201%2C%201%5D)%20//%209-,GPU%20%E2%80%94%203%20ms,-%40numba.cuda.jit @njit("void(double[:], double[:], double[:], double)") def bartlett_estimator_njit(result, k_norm, norm_xi_xj, order): for k in range(len(k_norm)): sum_k = 0.0 for ij in range(len(norm_xi_xj)): k_xi_xj = k_norm[k] * norm_xi_xj[ij] tmp = bessel1_njit(order, k_xi_xj) if order > 0: k_xi_xj = k_xi_xj ** order tmp /= k_xi_xj sum_k += tmp result[k] = 2 * sum_k