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)

sines = np.zeros_like(k)
widths_matrix = np.ones_like(k) * widths
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):
out = np.zeros_like(k)
widths_matrix = np.ones_like(k) * widths
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]