Source code for structure_factor.structure_factor

r"""Class collecting estimators of the structure factor :math:`S(\mathbf{k})` of stationary point process given one realization encapsulated in a :py:class:`~structure_factor.point_pattern.PointPattern`.

**The available estimators:**

- :py:meth:`~structure_factor.structure_factor.StructureFactor.scattering_intensity`: The scattering intensity and the corresponding debiased versions.

- :py:meth:`~structure_factor.structure_factor.StructureFactor.tapered_estimator`: The tapered or multitapered estimator and the corresponding debiased versions.

- :py:meth:`~structure_factor.structure_factor.StructureFactor.tapered_estimator_isotropic`: Bartlett's isotropic estimator.

- :py:meth:`~structure_factor.structure_factor.StructureFactor.quadrature_estimator_isotropic`: Integral estimation using Hankel transform quadrature.

**The available plot methods:**

- :py:meth:`~structure_factor.structure_factor.StructureFactor.plot_non_isotropic_estimator`: Visualize the results of :py:meth:`~structure_factor.structure_factor.StructureFactor.scattering_intensity`, :py:meth:`~structure_factor.structure_factor.StructureFactor.tapered_estimator`, or :py:meth:`~structure_factor.structure_factor.StructureFactor.multitapered_estimator`.

- :py:meth:`~structure_factor.structure_factor.StructureFactor.plot_isotropic_estimator`: Visualize the results of :py:meth:`~structure_factor.structure_factor.StructureFactor.tapered_estimator_isotropic` or :py:meth:`~structure_factor.structure_factor.StructureFactor.quadrature_estimator_isotropic`.

For the theoretical derivation and definitions of these estimators, we refer to :cite:`HGBLR:22`.
"""

import warnings

import numpy as np

import structure_factor.plotting as plots
from structure_factor.point_pattern import PointPattern
from structure_factor.spatial_windows import BallWindow, BoxWindow
from structure_factor.tapered_estimators import (
    allowed_k_scattering_intensity,
    scattering_intensity,
    select_tapered_estimator,
)
from structure_factor.tapered_estimators_isotropic import (
    allowed_k_norm_bartlett_isotropic,
    bartlett_estimator,
)
from structure_factor.tapers import BartlettTaper
from structure_factor.transforms import RadiallySymmetricFourierTransform


[docs]class StructureFactor: r"""Implementation of various estimators of the structure factor :math:`S(\mathbf{k})` of a stationary point process given one realization encapsulated in the :py:attr:`~structure_factor.structure_factor.StructureFactor.point_pattern` attribute. .. todo:: list attributes .. proof:definition:: The structure factor :math:`S` of a d-dimensional stationary point process :math:`\mathcal{X}` with intensity :math:`\rho` is defined by, .. math:: S(\mathbf{k}) = 1 + \rho \mathcal{F}(g-1)(\mathbf{k}), where :math:`\mathcal{F}` denotes the Fourier transform, :math:`g` the pair correlation function of :math:`\mathcal{X}`, :math:`\mathbf{k}` a wavevector of :math:`\mathbb{R}^d`. For more details we refer to :cite:`HGBLR:22`, (Section 2) or :cite:`Tor18`, (Section 2.1, equation (13)). """
[docs] def __init__(self, point_pattern): r"""Initialize StructureFactor from ``point_pattern``. Args: point_pattern (:py:class:`~structure_factor.point_pattern.PointPattern`): Object of type point pattern which contains a realization ``point_pattern.points`` of a point process, the window where the points were simulated ``point_pattern.window`` and (optionally) the intensity of the point process ``point_pattern.intensity``. """ assert isinstance(point_pattern, PointPattern) self.point_pattern = point_pattern
@property def dimension(self): """Ambient dimension of the underlying point process.""" return self.point_pattern.dimension
[docs] def scattering_intensity(self, k=None, debiased=True, direct=True, **params): r"""Compute the scattering intensity :math:`\widehat{S}_{\mathrm{SI}}` estimate of the structure factor from one realization of a stationary point process encapsulated in the :py:attr:`~structure_factor.structure_factor.StructureFactor.point_pattern` attribute. Args: k (numpy.ndarray, optional): Array of size :math:`n \times d` where :math:`d` is the dimension of the space, and :math:`n` is the number of wavevectors where the scattering intensity is evaluated. If ``k=None`` and ``debiased=True``, the scattering intensity will be evaluated on the corresponding set of allowed wavevectors; In this case, the keyword arguments ``k_max``, and ``meshgrid_shape`` can be used. Defaults to None. debiased (bool, optional): Trigger the use of a debiased tapered estimator. Defaults to True. If ``debiased=True``, the estimator is debiased as follows, - if ``k=None``, the scattering intensity will be evaluated on the corresponding set of allowed wavevectors. - if ``k`` is not None and ``direct=True``, the direct debiased scattering intensity will be used, - if ``k`` is not None and ``direct=False``, the undirect debiased scattering intensity will be used. direct (bool, optional): If ``debiased`` is True, trigger the use of the direct/undirect debiased scattering intensity. Parameter related to ``debiased``. Defaults to True. Keyword Args: params (dict): Keyword arguments ``k_max`` and ``meshgrid_shape`` of :py:func:`~structure_factor.tapered_estimators.allowed_k_scattering_intensity`, used when ``k=None`` and ``debiased=True``. Returns: tuple(numpy.ndarray, numpy.ndarray): - k: Wavevector(s) on which the scattering intensity has been evaluated. - estimation: Evaluation of the scattering intensity estimator of the structure factor at ``k``. Example: .. plot:: code/structure_factor/scattering_intensity.py :include-source: True .. proof:definition:: The scattering intensity :math:`\widehat{S}_{\mathrm{SI}}` is an estimator of the structure factor :math:`S` of a stationary point process :math:`\mathcal{X} \subset \mathbb{R}^d` with intensity :math:`\rho`. It is computed from one realization :math:`\mathcal{X}\cap W =\{\mathbf{x}_i\}_{i=1}^N` of :math:`\mathcal{X}` within a box window :math:`W=\prod_{j=1}^d[-L_j/2, L_j/2]`. .. math:: \widehat{S}_{\mathrm{SI}}(\mathbf{k}) = \frac{1}{N}\left\lvert \sum_{j=1}^N \exp(- i \left\langle \mathbf{k}, \mathbf{x_j} \right\rangle) \right\rvert^2 . For more details we refer to :cite:`HGBLR:22`, (Section 3.1). .. note:: **Typical usage** - If the observation window is not a :py:class:`~structure_factor.spatial_windows.BoxWindow`, use the method :py:meth:`~structure_factor.point_pattern.PointPattern.restrict_to_window` to extract a sub-sample in a :py:class:`~structure_factor.spatial_windows.BoxWindow`. .. seealso:: - :py:func:`~structure_factor.tapered_estimators.scattering_intensity` - :py:func:`~structure_factor.tapered_estimators.allowed_k_scattering_intensity` - :py:meth:`~structure_factor.structure_factor.StructureFactor.tapered_estimator` - :py:meth:`~structure_factor.structure_factor.StructureFactor.plot_non_isotropic_estimator` - :py:class:`~structure_factor.spatial_windows.BoxWindow` - :py:meth:`~structure_factor.point_pattern.PointPattern.restrict_to_window` """ if not isinstance(self.point_pattern.window, BoxWindow): warnings.warn( message="The observation window should be a BoxWindow to minimize the approximation error. Hint: use point_pattern.restrict_to_window." ) window = self.point_pattern.window d = window.dimension if k is None: if not debiased: raise ValueError("debiased argument must be True when k is None .") L = np.diff(window.bounds) k = allowed_k_scattering_intensity(d, L, **params) estimated_sf_k = scattering_intensity( k, self.point_pattern, debiased=debiased, direct=direct ) return k, estimated_sf_k
# todo consider unpacking *tapers to cover single/multiple taper/s at once?
[docs] def tapered_estimator(self, k, tapers, debiased=True, direct=True): r"""Compute the (multi)tapered estimator parametrized by ``tapers`` of the structure factor of a stationary point process given a realization encapsulated in ``point_pattern``. Args: k (numpy.ndarray): Array of size :math:`n \times d` where :math:`d` is the dimension of the space, and :math:`n` is the number of wavevectors where the tapered estimator is evaluated. t1, t2, ... : sequence of concrete tapers with methods ``.taper(x, window)`` corresponding to the taper function :math:`t(x, W)` , and ``.ft_taper(k, window)`` corresponding to the Fourier transform :math:`\mathcal{F}[t(\cdot, W)](k)` of the taper. See also :ref:`tapers`. debiased (bool, optional): Trigger the use of a debiased estimator. Defaults to True. direct (bool, optional): If ``debiased`` is True, trigger the use of the direct/undirect debiased tapered estimator. Parameter related to ``debiased``. Defaults to True. Returns: tuple(numpy.ndarray, numpy.ndarray): - k: Wavevector(s) on which the tapered estimator has been evaluated. - estimation: Evaluation of the tapered estimator at ``k``.. .. note:: Calling this method with its default arguments is equivalent to calling :py:meth:`~structure_factor.structure_factor.StructureFactor.scattering_intensity`. Example: .. plot:: code/structure_factor/multitapered_estimator.py :include-source: True .. proof:definition:: The tapered estimator :math:`\widehat{S}_{\mathrm{TP}}(t, k)`, is an estimator of the structure factor :math:`S` of a stationary point process :math:`\mathcal{X} \subset \mathbb{R}^d` with intensity :math:`\rho`. It is computed from one realization :math:`\mathcal{X}\cap W =\{\mathbf{x}_i\}_{i=1}^N` of :math:`\mathcal{X}` within a box window :math:`W`. .. math:: \widehat{S}_{\mathrm{TP}}(t, \mathbf{k}) = \frac{1}{\rho} \left\lvert \sum_{j=1}^N t(x_j, W) \exp(- i \left\langle k, x_j \right\rangle)\right\rvert^2, If several tapers are used, a simple average of the corresponding tapered estimators is computed. The resulting estimator is call multitapered estimator. .. math:: \widehat{S}_{\mathrm{MTP}}((t_{q})_{q=1}^P, \mathbf{k}) = \frac{1}{P}\sum_{q=1}^{P} \widehat{S}(t_{q}, \mathbf{k}) where, :math:`(t_{q})_{q}` is a family of tapers supported on the observation window (satisfying some conditions), :math:`P` is the number of tapers used, and :math:`k \in \mathbb{R}^d`. For more details, we refer to :cite:`HGBLR:22`, (Section 3.1). .. note:: **Typical usage** - If the observation window is not a :py:class:`~structure_factor.spatial_windows.BoxWindow`, use the method :py:meth:`~structure_factor.point_pattern.PointPattern.restrict_to_window` to extract a sub-sample in a :py:class:`~structure_factor.spatial_windows.BoxWindow`. .. seealso:: - :py:meth:`~structure_factor.structure_factor.StructureFactor.scattering_intensity` - :py:meth:`~structure_factor.structure_factor.StructureFactor.plot_non_isotropic_estimator` - :py:class:`~structure_factor.spatial_windows.BoxWindow` - :py:meth:`~structure_factor.point_pattern.PointPattern.restrict_to_window` - :ref:`tapers` - :py:class:`~structure_factor.tapers.SineTaper` - :py:func:`~structure_factor.tapered_estimators.multitapered_estimator` """ point_pattern = self.point_pattern n, d = k.shape if d != point_pattern.dimension: raise ValueError( f"k must be of size (n, d) where d=point_pattern.dimension. Given k {k.shape} and d d = {point_pattern.dimension}" ) _tapers = list(tapers) estimator = select_tapered_estimator(debiased, direct) estimated_sf_k = np.zeros(n, dtype=float) for t in _tapers: estimated_sf_k += estimator(k, point_pattern, t) estimated_sf_k /= len(_tapers) return k, estimated_sf_k
[docs] def bartlett_isotropic_estimator(self, k_norm=None, **params): r"""Compute Bartlett's isotropic estimator :math:`\widehat{S}_{\mathrm{BI}}` from one realization of an isotropic point process encapsulated in the :py:attr:`~structure_factor.structure_factor.StructureFactor.point_pattern` attribute. Args: k_norm (numpy.ndarray, optional): Array of wavenumbers of size :math:`n` where the estimator is to be evaluated. Defaults to None. Keyword Args: params (dict): Keyword argument ``nb_values`` of :py:func:`~structure_factor.tapered_estimators_isotropic.allowed_k_norm_bartlett_isotropic used when ``k_norm=None`` to specify the number of allowed wavenumbers to be considered. Returns: tuple(numpy.ndarray, numpy.ndarray): - k_norm: Wavenumber(s) on which Bartlett's isotropic estimator has been evaluated. - estimation: Evaluation(s) of Bartlett's isotropic estimator at ``k``. Example: .. plot:: code/structure_factor/bartlett_isotropic_estimator.py :include-source: True .. proof:definition:: Bartlett's isotropic estimator :math:`\widehat{S}_{\mathrm{BI}}` is an estimator of the structure factor :math:`S` of a stationary isotropic point process :math:`\mathcal{X} \subset \mathbb{R}^d` with intensity :math:`\rho`. It is computed from one realization :math:`\mathcal{X}\cap W =\{\mathbf{x}_i\}_{i=1}^N` of :math:`\mathcal{X}` within 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). .. note:: **Typical usage** - If the observation window is not a :py:class:`~structure_factor.spatial_windows.BallWindow`, use the method :py:meth:`~structure_factor.point_pattern.PointPattern.restrict_to_window` to extract a sub-sample in a :py:class:`~structure_factor.spatial_windows.BallWindow`. .. seealso:: - :py:class:`~structure_factor.spatial_windows.BallWindow` - :py:meth:`~structure_factor.point_pattern.PointPattern.restrict_to_window` - :py:func:`~structure_factor.tapered_estimators_isotropic` """ warnings.warn( message="The computation may take some time for a big number of points in the PointPattern. The complexity is quadratic in the number of points. Start by restricting the PointPattern to a smaller window using PointPattern.restrict_to_window, then increasing the window progressively." ) if not isinstance(self.point_pattern.window, BallWindow): warnings.warn( message="The observation window should be a BallWindow to minimize the approximation error. Hint: use point_pattern.restrict_to_window." ) if k_norm is None: if not isinstance(self.point_pattern.window, BallWindow): raise TypeError( "The observation window must be an instance of BallWindow. Hint: use point_pattern.restrict_to_window." ) window = self.point_pattern.window d, r = window.dimension, window.radius k_norm = allowed_k_norm_bartlett_isotropic( dimension=d, radius=r, **params ).astype(float) sf_k_norm = bartlett_estimator(k_norm, self.point_pattern) return k_norm, sf_k_norm
[docs] def quadrature_estimator_isotropic( self, pcf, k_norm=None, method="BaddourChouinard", **params ): # ? mettre k_nom avant pcf et donner le choix à l'utilisateur d'enter un None r"""Approximate the structure factor of a stationary isotropic point process at values ``k_norm``, given its pair correlation function ``pcf``, using a quadrature ``method``. .. warning:: It only applies to point processes in even dimension :math:`d`, due to evaluations of the zeros of Bessel functions of order :math:`d / 2 - 1` that must integer-valued. Args: pcf (callable): Pair correlation function. k_norm (numpy.ndarray, optional): Vector of wavenumbers (i.e., norms of wavevectors) where the structure factor is to be evaluated. Optional if ``method="BaddourChouinard"`` (since this method evaluates the Hankel transform on a specific set of wavenumbers), but it is **non optional** if ``method="Ogata"``. Defaults to None. method (str, optional): Trigger the use of ``"BaddourChouinard"`` or ``"Ogata"`` quadrature to estimate the structure factor. Defaults to ``"BaddourChouinard"``, - if ``"BaddourChouinard"``: The Hankel transform is approximated using the Discrete Hankel transform :cite:`BaCh15`. See :py:class:`~structure_factor.transforms.HankelTransformBaddourChouinard` - if ``"Ogata"``: The Hankel transform is approximated using Ogata quadrature :cite:`Oga05`. See :py:class:`~structure_factor.transforms.HankelTransformOgata` Keyword Args: params (dict): Keyword arguments passed to the corresponding Hankel transformer selected according to the input argument ``method``. - ``method == "Ogata"``, see :py:meth:`~structure_factor.transforms.HankelTransformOgata.compute_transformation_parameters` - ``step_size`` - ``nb_points`` - ``method == "BaddourChouinard"``, see :py:meth:`~structure_factor.transforms.HankelTransformBaddourChouinard.compute_transformation_parameters` - ``r_max`` - ``nb_points`` - ``interpolotation`` dictionnary containing the keyword arguments of `scipy.integrate.interp1d <https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html>`_ parameters. Returns: tuple (numpy.ndarray, numpy.ndarray): - k_norm: Vector of wavenumbers. - estimation: Evaluations of the structure factor on ``k_norm``. Example: .. plot:: code/structure_factor/hankel_quadrature.py :include-source: True .. proof:definition:: The structure factor :math:`S` of a **stationary isotropic** point process :math:`\mathcal{X} \subset \mathbb{R}^d` with intensity :math:`\rho`, can be defined via the Hankel transform :math:`\mathcal{H}_{d/2 -1}` of order :math:`d/2 -1` as follows, .. math:: S(\|\mathbf{k}\|_2) = 1 + \rho \frac{(2 \pi)^{d/2}}{\|\mathbf{k}\|_2^{d/2 -1}} \mathcal{H}_{d/2 -1}(\tilde g -1)(\|\mathbf{k}\|_2), \quad \tilde g: x \mapsto g(x) x^{d/2 -1}, where, :math:`g` is the pair correlation function of :math:`\mathcal{X}`. This is a result of the relation between the Symmetric Fourier transform and the Hankel Transform. For more details, we refer to :cite:`HGBLR:22`, (Section 3.2). .. note:: **Typical usage** 1. Estimate the pair correlation function using :py:meth:`~structure_factor.pair_correlation_function.PairCorrelationFonction.estimate`. 2. Clean and interpolate/extrapolate the resulting estimation using :py:meth:`~structure_factor.pair_correlation_function.PairCorrelationFonction.interpolate` to get a **function**. 3. Use the result as the input ``pcf``. .. seealso:: - :py:meth:`~structure_factor.pair_correlation_function.PairCorrelationFonction.estimate` - :py:meth:`~structure_factor.pair_correlation_function.PairCorrelationFonction.interpolate` - :py:meth:`~structure_factor.structure_factor.StructureFactor.plot_isotropic_estimator` - :py:class:`~structure_factor.spatial_windows` - :py:meth:`~structure_factor.point_pattern.PointPattern` - :py:class:`~structure_factor.transforms.HankelTransformBaddourChouinard` - :py:class:`~structure_factor.transforms.HankelTransformOgata` """ assert callable(pcf) if method == "Ogata" and k_norm.all() is None: raise ValueError( "k_norm argument must be passed when using method='Ogata'." ) r_max = params.setdefault("r_max", None) if method == "BaddourChouinard" and r_max is None: raise ValueError( "r_max keyword argument must be passed when using method='BaddourChouinard'." ) ft = RadiallySymmetricFourierTransform(dimension=self.dimension) total_pcf = lambda r: pcf(r) - 1.0 k_norm, ft_k = ft.transform(total_pcf, k_norm, method=method, **params) rho = self.point_pattern.intensity sf = 1.0 + rho * ft_k return k_norm, sf
[docs] def plot_isotropic_estimator( self, k_norm, estimation, axis=None, scale="log", k_norm_min=None, exact_sf=None, color="grey", error_bar=False, label=r"$\widehat{S}$", file_name="", **binning_params, ): r"""Display the outputs of the method :py:meth:`~structure_factor.structure_factor.StructureFactor.quadrature_estimator_isotropic`, or :py:meth:`~structure_factor.structure_factor.StructureFactor.tapered_estimator_isotropic`. Args: k_norm (numpy.ndarray): Vector of wavenumbers (i.e., norms of wavevectors) on which the structure factor has been approximated. estimation (numpy.ndarray): Approximation(s) of the structure factor corresponding to ``k_norm``. axis (plt.Axes, optional): Support axis of the plot. Defaults to None. scale (str, optional): Trigger between plot scales of `see matplolib documentation <https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.set_xscale.html>`_. Defaults to 'log'. k_norm_min (float, optional): Estimated lower bound of the wavenumbers. Defaults to None. exact_sf (callable, optional): Theoretical structure factor of the point process. Defaults to None. error_bar (bool, optional): If ``True``, ``k_norm`` and correspondingly ``estimation``, are divided into sub-intervals (bins). Over each bin, the mean and the standard deviation of ``si`` are derived and visualized on the plot. Note that each error bar corresponds to the mean :math:`\pm 3 \times` standard deviation. To specify the number of bins, add it as a keyword argument. For more details see :py:meth:`~structure_factor.utils._bin_statistics`. Defaults to False. file_name (str, optional): Name used to save the figure. The available output formats depend on the backend being used. Defaults to "". Keyword Args: binning_params: (dict): Used when ``error_bar=True``, by the method :py:meth:`~structure_factor.utils_bin_statistics` as keyword arguments (except ``"statistic"``) of `scipy.stats.binned_statistic <https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binned_statistic.html>`_. Returns: plt.Axes: Plot of the approximated structure factor. """ # todo normalize names of method and plotting routine return plots.plot_sf_hankel_quadrature( k_norm=k_norm, estimation=estimation, axis=axis, scale=scale, k_norm_min=k_norm_min, color=color, exact_sf=exact_sf, error_bar=error_bar, label=label, file_name=file_name, **binning_params, )
[docs] def plot_non_isotropic_estimator( self, k, estimation, axes=None, scale="log", plot_type="radial", positive=False, exact_sf=None, error_bar=False, label=r"$\widehat{S}$", rasterized=True, file_name="", window_res=None, **binning_params, ): r"""Display the outputs of the method :py:meth:`~structure_factor.structure_factor.StructureFactor.scattering_intensity`, :py:meth:`~structure_factor.structure_factor.StructureFactor.tapered_estimator`. Args: k (numpy.ndarray): Wavevector(s) on which the scattering intensity has been approximated. Array of size :math:`n \times d` where :math:`d` is the dimension of the space, and :math:`n` is the number of wavevectors. estimation (numpy.ndarray): Approximated structure factor associated to `k`. axes (plt.Axes, optional): Support axes of the plots. Defaults to None. scale (str, optional): Trigger between plot scales of `see matplolib documentation <https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.set_xscale.html>`_. Defaults to "log". plot_type (str, optional): Type of the plot to visualize, "radial", "imshow", or "all". Defaults to "radial". - If "radial", the output is a 1D plot of `estimation` w.r.t. the norm(s) of `k`. - If "imshow" (option available only for a 2D point process), the output is a 2D color level plot. - If "all" (option available only for a 2D point process), the result contains 3 subplots: the point pattern (or a restriction to a specific window if ``window_res`` is set), the radial plot, and the color level plot. Note that the options "imshow" and "all" couldn't be used, if ``k`` couldn't be reshaped as a meshgrid. positive (bool, optional): If True, consider only the positive values of `estimation`. Defaults to False. exact_sf (callable, optional): Theoretical structure factor of the point process. Defaults to None. error_bar (bool, optional): If ``True``, ``k_norm`` and correspondingly ``estimation``, are divided into sub-intervals (bins). Over each bin, the mean and the standard deviation of ``estimation`` are derived and visualized on the plot. Note that each error bar corresponds to the mean :math:`\pm 3 \times` standard deviation. To specify the number of bins, add it as a keyword argument. For more details see :py:meth:`~structure_factor.utils._bin_statistics`. Defaults to False. rasterized (bool, optional): Rasterized option of `matlplotlib.plot <https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html#:~:text=float-,rasterized,-bool>`_. Defaults to True. file_name (str, optional): Name used to save the figure. The available output formats depend on the backend being used. Defaults to "". window_res (:py:class:`~structure_factor.spatial_windows.AbstractSpatialWindow`, optional): New restriction window. Useful when the sample of points is large and "plt_type='all'", so for time and visualization purposes, it is better to restrict the plot of the point process to a smaller window. Defaults to None. Keyword Args: binning_params (dict): Used when ``error_bar=True``, by the method :py:meth:`~structure_factor.utils._bin_statistics` as keyword arguments (except ``"statistic"``) of `scipy.stats.binned_statistic <https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binned_statistic.html>`_. Returns: plt.Axes: Plot of the approximated structure factor. """ k_norm = np.linalg.norm(k, axis=1) # unplot possible negative values if positive: si_ = estimation estimation = estimation[si_ >= 0] k_norm = k_norm[si_ >= 0] if plot_type == "radial": return plots.plot_estimation_showcase( k_norm, estimation, axis=axes, scale=scale, exact_sf=exact_sf, error_bar=error_bar, label=label, rasterized=rasterized, file_name=file_name, **binning_params, ) elif plot_type == "imshow": if self.dimension != 2: raise ValueError( "This plot option is available only for a 2D point process. Hint: Use `plot_type ='radial'`." ) n_grid = int(np.sqrt(k_norm.shape[0])) if k_norm.shape[0] != n_grid ** 2: raise ValueError( "Wavevectors couldn't be reshaped as meshgrid. Hint: the square root of the number of wavevectors should be an integer." ) grid_shape = (n_grid, n_grid) estimation = estimation.reshape(grid_shape) k_norm = k_norm.reshape(grid_shape) return plots.plot_estimation_imshow(k_norm, estimation, axes, file_name) elif plot_type == "all": n_grid = int(np.sqrt(k_norm.shape[0])) grid_shape = (n_grid, n_grid) if self.dimension != 2: raise ValueError( "The option 'all' is available only for 2D point processes." ) estimation = estimation.reshape(grid_shape) k_norm = k_norm.reshape(grid_shape) return plots.plot_estimation_all( self.point_pattern, k_norm, estimation, exact_sf=exact_sf, error_bar=error_bar, label=label, rasterized=rasterized, file_name=file_name, window_res=window_res, scale=scale, **binning_params, ) else: raise ValueError( "plot_type must be chosen among ('all', 'radial', 'imshow')." )