Source code for structure_factor.transforms

"""Collections of classes that allow to compute the Fourier transform of a radially symmetric function and Hankel transforms.

- :py:class:`~structure_factor.transforms.RadiallySymmetricFourierTransform`: Compute the Fourier transform of a radially symmetric function using the `correspondence with the Hankel transform <https://en.wikipedia.org/wiki/Hankel_transform#Fourier_transform_in_d_dimensions_(radially_symmetric_case)>`_

- :py:class:`~structure_factor.transforms.HankelTransformBaddourChouinard`: Compute the Hankel transform using Baddour and Chouinard discrete Hankel transform

- :py:class:`~structure_factor.transforms.HankelTransformOgata`: Compute the Hankel transform using Ogata quadrature

For more details, we refer to :cite:`HGBLR:22`.
"""

import numpy as np
from scipy import interpolate

from structure_factor.utils import bessel1, bessel1_zeros, bessel2


[docs]class RadiallySymmetricFourierTransform: r"""Compute the Fourier transform of a radially symmetric function using the `correspondence with the Hankel transform <https://en.wikipedia.org/wiki/Hankel_transform#Fourier_transform_in_d_dimensions_(radially_symmetric_case)>`_. .. todo:: list attributes """
[docs] def __init__(self, dimension): """Initialize the :math:`d`-dimensional Fourier transform. Args: dimension (int): Dimension of the ambient space. """ assert isinstance(dimension, int) # assert dimension % 2 == 0 # required to evaluate zeros of Bessel functions with order = d / 2 - 1 that must integer-valued # error will be raised when calling bessel_zeros self.d = dimension self.r_max = None
[docs] def transform(self, f, k, method, **params): r"""Evaluate the Fourier transform of the radially symmetric function :math:`f` at :math:`k` using the correspondence with the Hankel transform. Args: f (callable): Function to transform. k (scalar or numpy.ndarray): Point or vector of points where the Fourier transform is to be evaluated. method (str): Name of the method used to compute the underlying Hankel transform. - ``"Ogata"`` :py:meth:`~structure_factor.transforms.HankelTransformOgata.transform` - ``"BaddourChouinard"`` :py:meth:`~structure_factor.transforms.HankelTransformBaddourChouinard.transform` Keyword Args: params (dict): - If ``method="BaddourChouinard"`` (see :cite:`BaCh15`): - r_max (float): Threshold radius characterizing the space-limited feature of the function ``f``, i.e., :math:`f(r)=0` for r > r_max. - nb_points (int, optional): Number of quadrature nodes. Defaults to 300. - see also :py:meth:`~structure_factor.transforms.HankelTransformBaddouChouinard.transform` - If ``method="Ogata"`` (see :cite:`Oga05`): - r_max (float, optional): Maximum radius on which the input function :math:`f` to be Hankel transformed was evaluated before the interpolation. Parameter used to conclude a lower bound on :math:`k` on which :math:`f` to be Hankel transformed. Defaults to None. - step_size (float, optional): Step size of the discretization scheme. Defaults to 0.01. - nb_points (int, optional): Number of quadrature nodes. Defaults to 300. - see also :py:meth:`~structure_factor.transforms.HankelTransformOgata.transform` Returns: tuple (numpy.ndarray, numpy.ndarray): - k: Point(s) where the Fourier transform is to be evaluated. - F_k: Fourier transform of ``f`` at ``k``. .. proof:definition:: The Hankel transform :math:`\mathcal{H}_{\nu}` of order :math:`\nu` of :math:`f` is defined by .. math:: \mathcal{H}_{\nu -1}(f)(k) = \int_0^\infty f(r) J_{\nu}(kr)r \mathrm{d}k, where :math:`J_{\nu}` is the Bessel function of first kind. The :math:`d`-dimensional Fourier transform :math:`\mathcal{F}` of the radially symmetric function :math:`f` at :math:`k` could be defined using the Hankel transform of :math:`x \rightarrow x^{d/2 -1}f(x)` of order :math:`d/2 -1` as follows, .. math:: k^{d/2-1} \mathcal{F}[f](k) = (2 \pi)^{d/2} \int_{0}^{+\infty} r^{d/2-1} f(r) J_{d/2-1}(kr) r \mathrm{d}r = (2 \pi)^{d/2} \mathcal{H}_{d/2-1}[\cdot^{d/2-1} f(\cdot)](k). """ d = self.d order = d // 2 - 1 ht = self._get_hankel_transformer(order, method) interp_params = params.pop("interpolation", dict()) ht.compute_transformation_parameters(**params) g = lambda r: f(r) * r ** order k, F_k = ht.transform(g, k, **interp_params) F_k *= (2 * np.pi) ** (d / 2) if order != 0: # F_k /= k^(d/2-1) F_k /= k ** order return k, F_k
def _get_hankel_transformer(self, order, method): hankel_transformer = { "Ogata": HankelTransformOgata, "BaddourChouinard": HankelTransformBaddourChouinard, } hankel_transform = hankel_transformer[method] return hankel_transform(order)
[docs]class HankelTransform: r"""Compute the `Hankel transform <https://en.wikipedia.org/wiki/Hankel_transform>`_ of order :math:`\nu`. .. todo:: list attributes .. seealso:: - :py:class:`~structure_factor.transforms.RadiallySymmetricFourierTransform` - :py:class:`~structure_factor.transforms.HankelTransformBaddourChouinard` - :py:class:`~structure_factor.transforms.HankelTransformOgata` """
[docs] def __init__(self, order): """Initialize the Hankel transform with prescribed ``order``. Args: order (int, optional): Order of the Hankel transform. """ assert order == np.floor(order) self.order = int(order)
[docs]class HankelTransformBaddourChouinard(HankelTransform): r"""Compute the Hankel transform, using the method of :cite:`BaCh15` considering that the input function is space-limited, i.e., :math:`f(r)=0` for :math:`r>r_{max}`. .. todo:: list attributes .. seealso:: - `MatLab code of Baddour Chouinard <https://openresearchsoftware.metajnl.com/articles/10.5334/jors.82/>`_ - `Pyhank Python package <https://pypi.org/project/pyhank/>`_ """
[docs] def __init__(self, order=0): """Initialize the Hankel transform with prescribed ``order``. Args: order (int, optional): Order of the Hankel transform. """ super().__init__(order=order) self.bessel_zeros = None self.r_max = None # R in :cite:`BaCh15` Section 4.B self.transformation_matrix = None # Y in :cite:`BaCh15` Section 6.A
[docs] def compute_transformation_parameters(self, r_max, nb_points): r"""Compute the parameters involved in the evaluation of the corresponding Hankel-type transform using the discretization scheme of :cite:`BaCh15`. The following object's attributes are defined - :py:attr:`~structure_factor.transforms.HankelTransformBaddourChouinard.bessel_zeros` - :py:attr:`~structure_factor.transforms.HankelTransformBaddourChouinard.r_max` - :py:attr:`~structure_factor.transforms.HankelTransformBaddourChouinard.transformation_matrix` Args: r_max (float): Threshold radius. Considering that the input function :math:`f` to be Hankel transformed is space-limited, then ``r_max`` satisfies :math:`f(r)=0` for r > r_max. nb_points (int): Number of quadrature nodes. """ n = self.order bessel_zeros = bessel1_zeros(n, nb_points) jk, jN = bessel_zeros[:-1], bessel_zeros[-1] # Section 6.A Transformation matrix Y = bessel1(n, np.outer(jk / jN, jk)) / np.square(bessel1(n + 1, jk)) Y *= 2 / jN self.bessel_zeros = bessel_zeros self.r_max = r_max self.transformation_matrix = Y
[docs] def transform(self, f, k=None, **interpolation_params): r"""Compute the Hankel transform of ``f`` at ``k``. Args: f (callable): Function to be Hankel transformed. k (numpy.ndarray, optional): Points of evaluation of the Hankel transform. Defaults to None. - If ``k`` is None (default), then ``k = self.bessel_zeros[:-1] / self.r_max`` derived from :py:meth:`~structure_factor.transforms.HankelTransformBaddourChouinard.compute_transformation_parameters`. - If ``k`` is provided, the Hankel transform is first computed at the above k values (case k is None), then interpolated using :py:func:`scipy.interpolate.interp1d` with ``interpolation_params`` and finally evaluated at the provided ``k`` values. Keyword Args: interpolation_params (dict): Keyword arguments of :py:func:`scipy.interpolate.interp1d`. Returns: tuple (scalar or numpy.ndarray, scalar or numpy.ndarray): ``k`` and the evaluations of the Hankel transform of ``f`` at ``k``. """ assert callable(f) r_max = self.r_max Y = self.transformation_matrix jk, jN = self.bessel_zeros[:-1], self.bessel_zeros[-1] r = jk * (r_max / jN) ht_k = (r_max ** 2 / jN) * Y.dot(f(r)) # Equation (23) _k = jk / r_max if k is None: return _k, ht_k interpolation_params["assume_sorted"] = True interpolation_params.setdefault("fill_value", "extrapolate") interpolation_params.setdefault("kind", "cubic") ht = interpolate.interp1d(_k, ht_k, **interpolation_params) return k, ht(k)
[docs]class HankelTransformOgata(HankelTransform): r"""Compute the Hankel transform using Ogata quadrature :cite:`Oga05`, (Section 5). .. todo:: list attributes .. seealso:: - `Hankel Python package <https://joss.theoj.org/papers/10.21105/joss.01397>`_ """
[docs] def __init__(self, order=0): """Initialize the Hankel transform with prescribed ``order``. Args: order (int, optional): Order of the Hankel transform. Defaults to 0. """ super().__init__(order=order) self.nodes, self.weights = None, None
[docs] def compute_transformation_parameters( self, r_max=None, nb_points=300, step_size=0.01 ): """Compute the quadrature nodes and weights used by :cite:`Oga05` (Equation (5.2)), to evaluate the corresponding Hankel-type transform. Args: r_max (float, optional): Maximum radius on which the input function :math:`f`, to be Hankel transformed, was evaluated before the interpolation. Parameter used to conclude a lower bound on :math:`k` on which :math:`f` to be Hankel transformed. Defaults to None. step_size (float, optional): Step size of the discretization scheme. Defaults to 0.01. nb_points (int, optional): Number of quadrature nodes. Defaults to 300. Returns: tuple (numpy.ndarray, np.ndarray): Quadrature nodes and weights. """ n = self.order h = step_size N = nb_points self.r_max = r_max t = bessel1_zeros(n, N) weights = bessel2(n, t) / bessel1(n + 1, t) # Equation (1.2) t *= h / np.pi # Equivalent of xi variable weights *= self._d_psi(t) nodes = (np.pi / h) * self._psi(t) # Change of variable Equation (5.1) self.nodes, self.weights = nodes, weights return nodes, weights
[docs] def transform(self, f, k): r"""Compute the Hankel transform of ``f`` evaluated at ``k``, following the work of :cite:`Oga05` (Section 5). Args: f (callable): Function to be Hankel transformed. k (numpy.ndarray, optional): Points of evaluation of the Hankel transform (1d array). Defaults to None. Returns: tuple(numpy.ndarray, np.ndarray): - k: Points of evaluation of the Hankel transform. - H_k: Evaluations of the Hankel transform of ``f`` on ``k``. .. important:: Please call :py:meth:`HankelTransformOgata.compute_transformation_parameters` to define quadrature attributes :py:attr:`~structure_factor.transforms.Ogata.nodes` and :py:attr:`~structure_factor.transforms.Ogata.weights`, before applying :py:meth:`HankelTransformOgata.compute_transform`. """ assert callable(f) n = self.order w = self.weights x = self.nodes k_ = k[:, None] if isinstance(k, np.ndarray) else k g = lambda r: f(r / k_) * r # or f(r / k_) * (r / k**2) H_k = np.pi * np.sum(w * g(x) * bessel1(n, x), axis=-1) H_k /= k ** 2 return k, H_k
@staticmethod def _psi(t): """Change of variable used by :cite:`Oga05` Equation (5.1).""" return t * np.tanh((0.5 * np.pi) * np.sinh(t)) @staticmethod def _d_psi(t): """Change of variable used by :cite:`Oga05` Equation (5.1).""" threshold = 3.5 # threshold outside of which psi' plateaus to -1, 1 out = np.sign(t) mask = np.abs(t) < threshold x = t[mask] out[mask] = np.pi * x * np.cosh(x) + np.sinh(np.pi * np.sinh(x)) out[mask] /= 1.0 + np.cosh(np.pi * np.sinh(x)) return out