Source code for structure_factor.tapers

r"""Collection of classes representing tapers (also called tapering functions or window functions). Such tapers must have two methods.

- ``.taper(x, window)`` corresponding to the tapering function :math:`t(x, W)`,

- ``.ft_taper(k, window)`` corresponding to the Fourier transform :math:`\mathcal{F}[t(\cdot, W)](k)` of the tapering function.

These tapers satisfy some conditions listed in :cite:`HGBLR:22` (Sections 3.1, 4.3).

Example:
    .. literalinclude:: code/structure_factor/taper.py
        :language: python
"""

from itertools import product

import numpy as np

from structure_factor.spatial_windows import BoxWindow


[docs]class BartlettTaper: """Class representing the Bartlett tapering function."""
[docs] @staticmethod def taper(x, window): r"""Evaluate the Bartlett taper :math:`t(x, W)` at ``x`` given the rectangular ``window`` :math:`W`. .. math:: t(x, W) = \frac{1}{\sqrt{|W|}} 1_{x \in W}. Args: x (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 tapering function is evaluated. window (:py:class:`~structure_factor.spatial_windows.BoxWindow`): :math:`d`-dimensional rectangular window :math:`W`. Returns: numpy.ndarray: evaluation of the taper :math:`t(x, W)`. """ assert isinstance(window, BoxWindow) t = window.indicator_function(x).astype(float) t /= np.sqrt(window.volume) return t
[docs] @staticmethod def ft_taper(k, window): r"""Evaluate the Fourier transform :math:`F[t(\cdot, W)](k)` of the taper :math:`t` (:py:meth:`~structure_factor.tapers.BartlettTaper.taper`). 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 Fourier transform is evaluated. window (:py:class:`~structure_factor.spatial_windows.BoxWindow`): :math:`d`-dimensional rectangular window :math:`W`. Return: numpy.ndarray: Evaluation of the Fourier transform at ``k``. """ assert isinstance(window, BoxWindow) widths = 0.5 * np.diff(window.bounds.T, axis=0) mask = k == 0 if np.any(mask): sines = np.zeros_like(k) widths_matrix = np.ones_like(k) * widths sines[mask] = 2.0 * widths_matrix[mask] np.logical_not(mask, out=mask) sines[mask] = 2.0 * np.sin(k[mask] * widths_matrix[mask]) / k[mask] else: sines = 2.0 * np.sin(k * widths) / k ft = np.prod(sines, axis=1) return ft / np.sqrt(window.volume)
[docs]class SineTaper: r"""Class representing the sine tapering function."""
[docs] def __init__(self, p): self.p = np.array(p)
[docs] def taper(self, x, window): r"""Evalute the sine taper :math:`t(x, W)` indexed by :py:attr:`~structure_factor.tapers.SineTaper.p` at ``x`` given the rectangular ``window`` :math:`W`. Args: x (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 tapering function is evaluated. window (:py:class:`~structure_factor.spatial_windows.BoxWindow`): :math:`d`-dimensional rectangular window :math:`W`. Returns: numpy.ndarray: evaluation of the taper :math:`t(x, W)`. """ assert isinstance(window, BoxWindow) widths = np.diff(window.bounds.T, axis=0) # d = x.shape[1] sines = x / widths + 0.5 sines *= np.pi * self.p np.sin(sines, out=sines) t = window.indicator_function(x).astype(float) t *= np.prod(sines * np.sqrt(2), axis=1) t *= np.sqrt(1 / window.volume) # normalization return t
[docs] def ft_taper(self, k, window): r"""Evaluate the Fourier transform :math:`F[t(\cdot, W)](k)` of the taper :math:`t` (:py:meth:`~structure_factor.tapers.SineTaper.taper`). 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 Fourier transform is evaluated. window (:py:class:`~structure_factor.spatial_windows.BoxWindow`): :math:`d`-dimensional rectangular window :math:`W`. Return: numpy.ndarray: Evaluation of the Fourier transform at ``k``. """ assert isinstance(window, BoxWindow) widths = np.diff(window.bounds.T, axis=0) p = self.p a = k - np.pi * p / widths b = k + np.pi * p / widths def compute_factor(x, k, widths): mask = x == 0 if np.any(mask): out = np.zeros_like(k) widths_matrix = np.ones_like(k) * widths out[mask] = 0.5 * widths_matrix[mask] np.logical_not(mask, out=mask) out[mask] = np.sin(x[mask] * 0.5 * widths_matrix[mask]) / x[mask] else: out = np.sin(x * 0.5 * widths) / x return out res_a = compute_factor(a, k, widths) res_b = compute_factor(b, k, widths) res = (1j) ** (p + 1) * np.sqrt(2) * (res_a - (-1) ** p * res_b) ft = np.prod(res, axis=1) ft /= np.sqrt(window.volume) return ft
[docs]def multi_sinetaper_grid(d, p_component_max=2): r"""Given a class of taper `taper_p` of parameter `p` of :math:`\mathbb{R}^d`, return the list of taper `taper_p(p)` with :math:`p \in \{1, ..., P\}^d`. Args: d (int): Space dimension. taper_p (Class): Class of taper pf parameter p. p_component_max (int): Maximum component of the parameters :math:`p` of the family of tapers. Intuitively the number of taper used is :math:`P=\mathrm{p\_component\_max}^d`. Used only when ``tapers=None``. Defaults to 2. Returns: list: List of taper `taper_p(p)` with :math:`p \in \{1, ..., p\_component\_max\}^d`. """ params = product(*(range(1, p_component_max + 1) for _ in range(d))) return [SineTaper(p) for p in params]