Source code for structure_factor.multiscale_estimators

import warnings
import numpy as np
from structure_factor.structure_factor import StructureFactor
from structure_factor.spatial_windows import (
from structure_factor.tapered_estimators_isotropic import (
from scipy.stats import poisson

[docs]def multiscale_estimator( point_pattern, estimator, subwindows_list, k_list, mean_poisson, m=None, proba_list=None, verbose=True, **kwargs, ): r"""Sample from :math:`Z` :cite:`HGBLR:22` using a PointPattern and a realization from the r.v. :math:`M`. See the definition of :math:`Z` below. Args: point_pattern (:py:class:`~structure_factor.point_pattern.PointPattern`): An encapsulation of a realization of a point process, the observation window, and (optionally) the intensity of the point process. 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 :py:class:`~structure_factor.structure_factor.StructureFactor`. subwindows_list (list): List of increasing cubic or ball-shaped :py:class:`~structure_factor.spatial_windows.AbstractSpatialWindow`, typically, obtained using :py:func:`~structure_factor.multiscale_estimators.subwindows_list`. 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 :math:`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 :py:func:`~structure_factor.multiscale_estimators.subwindows_list`. mean_poisson (int): Parameter of the Poisson r.v. :math:`M` used to compute :math:`Z`. To use a different distribution of the r.v. :math:`M`, set ``mean_poisson=None`` and specify ``m_list`` and ``proba_list`` corresponding to :math:`M`. m (int, optional): Realization of the positive integer-valued r.v. :math:`M` used when ``mean_poisson=None``. Defaults to None. proba_list (list, optional): List of :math:`\mathbb{P}(M \geq j)` used with ``m`` when ``mean_poisson=None``. Should contains at least ``m`` elements. Defaults to None. verbose (bool, optional): If "True" and ``mean_poisson`` is not None, print the re-sampled values of :math:`M`. Defaults to False. Keyword Args: kwargs (dict): Parameters of the chosen ``estimator`` of the structure factor. See :py:class:`~structure_factor.structure_factor.StructureFactor`. Returns: float: The obtained value of :math:`Z`. Example: .. literalinclude:: code/multiscale_estimators/ :language: python .. proof:definition:: Let :math:`\mathcal{X} \in \mathbb{R}^d` be a stationary point process of which we consider an increasing sequence of sets :math:`(\mathcal{X} \cap W_m)_{m \geq 1}`, with :math:`(W_m)_m` centered box (or ball)-shaped windows s.t. :math:`W_s \subset W_r` for all :math:`0< s<r`, and :math:`W_{\infty} = \mathbb{R}^d`. We define the sequence of r.v. :math:`Y_m = 1\wedge \widehat{S}_m(\mathbf{k}_m^{\text{min}})`, where :math:`\widehat{S}_m` is one of the positive, asymptotically unbiased estimators of the structure factor of :math:`\mathcal{X}` applied on the observation :math:`\mathcal{X} \cap W_m`, and :math:`\mathbf{k}_m^{\text{min}}` is the minimum allowed wavevector associated with :math:`W_m`. Then, :math:`Z` is defined by, .. math:: Z = \sum_{j=1}^{M} \frac{Y_j - Y_{j-1}}{\mathbb{P}(M\geq j)} with :math:`M` an :math:`\mathbb{N}`-valued random variable such that :math:`\mathbb{P}(M \geq j)>0` for all :math:`j`, and :math:`Y_{0}=0` :cite:`HGBLR:22`, :cite:`RhGl15`. """ # r.v. threshold m_thresh = m_threshold( window_min=subwindows_list[0], window_max=subwindows_list[-1], ) # r.v. M if m is None: m = _poisson_rv(mean_poisson, m_thresh, verbose=verbose) else: if m > m_thresh: warnings.warn( message=f"The random variable M exceed the allowed threshold {m_thresh}." ) m = int(m) # proba list if proba_list is None: proba_list = 1 - ( poisson.cdf(k=range(1, m + 1), mu=mean_poisson) - poisson.pmf(k=range(1, m + 1), mu=mean_poisson) ) else: if len(proba_list) < m: raise ValueError(f"The proba list should contains at least{m} elements.") proba_list = proba_list[:m] # k and subwindows list if len(subwindows_list) != len(k_list): raise ValueError( "The number of wavevectors/wavenumber (k) should be equal to the number of subwindows. Each k is associated to a subwindow." ) if len(subwindows_list) < m: raise ValueError( f"The number of subwindows {len(subwindows_list)} should be at least equal to the random variable M= {m}." ) subwindows_list = subwindows_list[:m] k_list = k_list[:m] # approximated s_k_min list s_k_min_list = multiscale_estimator_core( point_pattern=point_pattern, subwindows_list=subwindows_list, k_list=k_list, estimator=estimator, **kwargs, ) # the r.v. (y_n)_n y_list = [min(np.array([1]), s) for s in s_k_min_list] z = coupled_sum_estimator(y_list, proba_list) return z
def multiscale_estimator_core( point_pattern, subwindows_list, k_list, estimator, **kwargs ): point_pattern_list = [point_pattern.restrict_to_window(w) for w in subwindows_list] estimated_sf_k_list = [ _select_structure_factor_estimator( point_pattern=p, estimator=estimator, k=q, **kwargs ) for p, q in zip(point_pattern_list, k_list) ] return estimated_sf_k_list
[docs]def coupled_sum_estimator(y_list, proba_list): r"""The coupled sum estimator of :cite:`RhGl15`. Args: y_list (list): List of :math:`M` realizations of the r.v. :math:`Y`. proba_list (list): List of :math:`\mathbb{P}(M \geq j)` with :math:`1 \leq j \leq M`. Returns: float: Obtained value of the coupled sum estimator. Example: .. literalinclude:: code/multiscale_estimators/ :language: python .. proof:definition:: Let :math:`(Y_m)_{m\geq 1}` be a sequence of :math:`L^2` approximations of a r.v. :math:`Y` each of which can be generated in finite time, for which :math:`\mathbb{E}^{1/2}[(Y_m - Y)^2]` converges to zero as :math:`m` goes to infinity. The coupled sum estimator of :cite:`RhGl15` is defined by, .. math:: Z = \sum_{j=1}^{M} \frac{Y_j - Y_{j-1}}{\mathbb{P}(M\geq j)}, with :math:`M` an :math:`\mathbb{N}`-valued random variable such that :math:`\mathbb{P}(M \geq j)>0` for all :math:`j`, and :math:`Y_{0}=0`. """ y_list_with_0 = np.append(0, y_list) # 0 first element of the list y_pairwise_diff = np.array( [t - s for s, t in zip(y_list_with_0[:-1], y_list_with_0[1:])] ) y_pairwise_diff = y_pairwise_diff / np.array(proba_list) z = np.sum(y_pairwise_diff) return z
[docs]def subwindows_list( window, subwindows_type="BoxWindow", param_0=None, param_max=None, params=None ): """Create a list of cubic (or ball)-shaped subwindows of a father window, with the associated minimum allowed wavevectors (or wavenumbers). Args: window (:py:class:`~structure_factor.spatial_windows.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: (list, list): - subwindows: Obtained subwindows. - k: Minimum allowed wavevectors of :py:func:`~structure_factor.tapered_estimators.allowed_k_scattering_intensity` or wavenumbers of :py:func:`~structure_factor.tapered_estimators_isotropic.allowed_k_norm_bartlett_isotropic` associated with the subwindow list. The former is for cubic and the latter for ball-shaped subwindows. Example: .. plot:: code/multiscale_estimators/ :include-source: True .. seealso:: - :py:mod:`~structure_factor.spatial_windows` - :py:func:`~structure_factor.tapered_estimators.allowed_k_scattering_intensity` - :py:func:`~structure_factor.tapered_estimators_isotropic.allowed_k_norm_bartlett_isotropic` """ d = window.dimension # parameter of the biggest possible subwindow contained in `window`` window_param_max = subwindow_parameter_max(window, subwindows_type) # subwindows list of parameters if params is None: if param_0 is None: raise ValueError( "The minimum window parameter is mandatory. Hint: specify the minimum window parameter." ) # from param_0 till param_max with unit space if param_max is None: params = np.arange(param_0, window_param_max) else: if max(params) > window_param_max: raise ValueError( f"The maximum sub-window (parameter={max(params)}) is bigger than the father window (parameter={window_param_max}). Hint: Reduce the maximum subwindow parameter. " ) # subwindows and the associated k if subwindows_type == "BallWindow": subwindows = [BallWindow(center=[0] * d, radius=r) for r in params] # check if d is even _, mod_ = divmod(d, 2) is_even_dimension = mod_ == 0 if is_even_dimension: k_list = [ allowed_k_norm_bartlett_isotropic(dimension=d, radius=r, nb_values=1) for r in params ] else: k_list = None warnings.warn( message=f"Isotropic allowed wavenumber are available only when the dimension of the space is an even number (i.e., d/2 is an integer)." ) else: subwindows = [BoxWindow(bounds=[[-l / 2, l / 2]] * d) for l in params] k_list = [np.full((1, d), fill_value=2 * np.pi / l) for l in params] return subwindows, k_list
[docs]def m_threshold(window_min, window_max): r"""Find the maximum number of integers ranging between the parameters (lengthside/radius) of ``window_min`` and the largest subwindow of ``window_max`` having the shape of ``window_min``. In particular, it gives the maximum value of the r.v. :math:`M` that can be used to compute the :py:func:`~structure_factor.multiscale_estimators.multiscale_estimator` given the smallest and biggest subwindow. Args: window_min (:py:class:`~structure_factor.spatial_windows.AbstractSpatialWindow`): Smallest cubic or ball-shaped window centered at the origin. window_max (:py:class:`~structure_factor.spatial_windows.AbstractSpatialWindow`): Biggest box or ball-shaped window centered at the origin. Returns: int : Maximum number of integers ranging between the parameters (lengthside/radius) of ``window_min`` and the largest subwindow of ``window_max`` having the shape of ``window_min``. Example: .. plot:: code/multiscale_estimators/ :include-source: True .. seealso:: - :py:mod:`~structure_factor.spatial_windows` """ if isinstance(window_min, BoxWindow): subwindow_type = "BoxWindow" else: subwindow_type = "BallWindow" param_max = subwindow_parameter_max(window_max, subwindow_type=subwindow_type) param_min = subwindow_parameter_max(window_min, subwindow_type=subwindow_type) if param_max < param_min: raise ValueError("window_min should be bigger than window_max.") m_threshold = int(param_max - param_min) return m_threshold
def _select_structure_factor_estimator(point_pattern, estimator, k, **kwargs): sf = StructureFactor(point_pattern) if estimator == "scattering_intensity": _, estimated_sf_k = sf.scattering_intensity(k=k, **kwargs) elif estimator == "tapered_estimator": _, estimated_sf_k = sf.tapered_estimator(k=k, **kwargs) elif estimator == "bartlett_isotropic_estimator": _, estimated_sf_k = sf.bartlett_isotropic_estimator(k_norm=k, **kwargs) elif estimator == "quadrature_estimator_isotropic": _, estimated_sf_k = sf.quadrature_estimator_isotropic(k_norm=k, **kwargs) else: raise ValueError( "Available estimators: 'scattering_intensity', 'tapered_estimator', 'bartlett_isotropic_estimator', 'quadrature_estimator_isotropic'. " ) return estimated_sf_k def _poisson_rv(mean_poisson, threshold, verbose=True): m = poisson.rvs(mu=mean_poisson, size=1) while m > threshold: if verbose: print("Re-sample M; current M= ", m, ", threshold=", threshold) m = int(poisson.rvs(mu=mean_poisson, size=1)) return m