diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..8ba354a --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +__pycache__ +*.swp +*.egg-info +dist + +docs/reference +docs/_build diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..6348618 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,5 @@ +repos: +- repo: https://github.com/pycqa/flake8 + rev: 4.0.1 + hooks: + - id: flake8 \ No newline at end of file diff --git a/.readthedocs.yml b/.readthedocs.yml new file mode 100644 index 0000000..56bcdef --- /dev/null +++ b/.readthedocs.yml @@ -0,0 +1,8 @@ +version: 2 + +python: + install: + - method: pip + path: . + extra_requirements: + - docs \ No newline at end of file diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..c82970e --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2022 Nicolas Tessore + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 index 0000000..d3fbd68 --- /dev/null +++ b/README.md @@ -0,0 +1,15 @@ +`fftl` — generalised FFTLog for Python +====================================== + +The `fftl` package for Python contains a routine to calculate integral +transforms of the type *ã(k) = ∫ a(r) T(kr) dr* for arbitrary kernels *T*. It +uses a modified FFTLog[^2] method of Hamilton[^1] to efficiently compute the +transform on logarithmic input and output grids. + +The package only requires `numpy`. To install with `pip`: + + pip install fftl + + +[^1]: Hamilton A. J. S., 2000, MNRAS, 312, 257 (astro-ph/9905191) +[^2]: Talman J. D., 1978, J. Comp. Phys., 29, 35 diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d4bb2cb --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..1f836a2 --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,68 @@ +# Configuration file for the Sphinx documentation builder. +# +# This file only contains a selection of the most common options. For a full +# list see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys +sys.path.insert(0, os.path.abspath('..')) + + +# -- Project information ----------------------------------------------------- + +project = 'fftl' +copyright = '2022, Nicolas Tessore' +author = 'Nicolas Tessore' + + +# -- General configuration --------------------------------------------------- + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'sphinx.ext.autosummary', + 'numpydoc', + 'matplotlib.sphinxext.plot_directive', +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'furo' + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = [] + + +# -- numpydoc extension ----------------------------------------------------- + +numpydoc_use_plots = True +numpydoc_show_class_members = False + + +# -- matplotlib extension ---------------------------------------------------- + +plot_include_source = True +plot_html_show_source_link = False +plot_html_show_formats = False diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 index 0000000..2459549 --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,37 @@ +***************************************** +``fftl`` -- generalised FFTLog for Python +***************************************** + +.. currentmodule:: fftl + +The ``fftl`` package for Python contains a routine to calculate integral +transforms of the type + +.. math:: + + \tilde{a}(k) = \int_{0}^{\infty} \! a(r) \, T(kr) \, dr + +for arbitrary kernels :math:`T`. It uses a modified FFTLog [2]_ method of +Hamilton [1]_ to efficiently compute the transform on logarithmic input and +output grids. + +The package only requires ``numpy``. To install with ``pip``:: + + pip install fftl + + +References +========== + +.. [1] Hamilton A. J. S., 2000, MNRAS, 312, 257 (astro-ph/9905191) +.. [2] Talman J. D., 1978, J. Comp. Phys., 29, 35 + + +Functions +========= + +.. autosummary:: + :toctree: reference + :nosignatures: + + fftl diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000..32bb245 --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=. +set BUILDDIR=_build + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/fftl.py b/fftl.py new file mode 100644 index 0000000..fd9b5b5 --- /dev/null +++ b/fftl.py @@ -0,0 +1,223 @@ +# author: Nicolas Tessore +# license: MIT +'''generalised FFTLog module''' + +__version__ = '2022.7.13' + +import numpy as np + +PI = np.pi + + +def fftl(u, r, ar, q=0.0, args=(), kr=1.0, krgood=True, deriv=False): + r'''generalised FFTLog for integral transforms + + Computes integral transforms for arbitrary kernels using a generalisation + of Hamilton's method [1]_ for the FFTLog algorithm [2]_. + + The kernel of the integral transform is characterised by the coefficient + function ``u``, see notes below, which must be callable and accept complex + input arrays. Additional arguments for ``u`` can be passed with the + optional ``args`` parameter. + + The function to be transformed must be given on a logarithmic grid ``r``. + The result of the integral transform is similarly computed on a logarithmic + grid ``k = kr/r``, where ``kr`` is a scalar constant (default: 1) which + shifts the logarithmic output grid. The selected value of ``kr`` is + automatically changed to the nearest low-ringing value if ``krgood`` is + true (the default). + + The integral transform can optionally be biased, see notes below. + + The function can optionally at the same time return the derivative of the + integral transform with respect to the logarithm of ``k``, by setting + ``deriv`` to true. + + Parameters + ---------- + u : callable + Coefficient function. Must have signature ``u(x, *args)`` and support + complex input arrays. + r : array_like (N,) + Grid of input points. Must have logarithmic spacing. + ar : array_like (..., N) + Function values. If multidimensional, the integral transform applies + to the last axis, which must agree with input grid. + q : float, optional + Bias parameter for integral transform. + args : tuple, optional + Additional arguments for the coefficient function ``u``. + kr : float, optional + Shift parameter for logarithmic output grid. + krgood : bool, optional + Change given ``kr`` to the nearest value fulfilling the low-ringing + condition. + deriv : bool, optional + Also return the first derivative of the integral transform. + + Returns + ------- + k : array_like (N,) + Grid of output points. + ak : array_like (..., N) + Integral transform evaluated at ``k``. + dak : array_like (..., N), optional + If ``deriv`` is true, the derivative of ``ak`` with respect to the + logarithm of ``k``. + + Notes + ----- + Computes integral transforms of the form + + .. math:: + + \tilde{a}(k) = \int_{0}^{\infty} \! a(r) \, T(kr) \, dr + + for arbitrary kernels :math:`T`. + + If :math:`a(r)` is given on a logarithmic grid of :math:`r` values, the + integral transform can be computed for a logarithmic grid of :math:`k` + values with a modification of Hamilton's FFTLog algorithm, + + .. math:: + + U(x) = \int_{0}^{\infty} \! t^x \, T(t) \, dt \;. + + The generalised FFTLog algorithm therefore only requires the coefficient + function :math:`U` for the given kernel. Everything else, and in + particular how to construct a well-defined transform, remains exactly the + same as in Hamilton's original algorithm. + + The transform can optionally be biased, + + .. math:: + + \tilde{a}(k) = k^{-q} \int_{0}^{\infty} \! [a(r) \, r^{-q}] \, + [T(kr) \, (kr)^q] \, dr \;, + + where :math:`q` is the bias parameter. The respective biasing factors + :math:`r^{-q}` and :math:`k^{-q}` for the input and output values are + applied internally. + + References + ---------- + .. [1] Hamilton A. J. S., 2000, MNRAS, 312, 257 (astro-ph/9905191) + .. [2] Talman J. D., 1978, J. Comp. Phys., 29, 35 + + Examples + -------- + Compute the one-sided Laplace transform of the hyperbolic tangent function. + The kernel of the Laplace transform is :math:`\exp(-kt)`, which determines + the coefficient function. + + >>> from scipy.special import gamma, digamma + >>> + >>> def u_laplace(x): + ... # requires Re(x) = q > -1 + ... return gamma(1 + x) + + Create the input function values on a logarithmic grid. + + >>> r = np.logspace(-4, 4, 100) + >>> ar = np.tanh(r) + >>> + >>> import matplotlib.pyplot as plt + >>> plt.loglog(r, ar) + >>> plt.xlabel('$r$') + >>> plt.ylabel('$\\tanh(r)$') + >>> plt.show() + + Compute the Laplace transform, and compare with the analytical result. + + >>> from fftl import fftl + >>> + >>> k, ak = fftl(u_laplace, r, ar) + >>> + >>> lt = (digamma((k+2)/4) - digamma(k/4) - 2/k)/2 + >>> + >>> plt.loglog(k, ak) + >>> plt.loglog(k, lt, ':') + >>> plt.xlabel('$k$') + >>> plt.ylabel('$L[\\tanh](k)$') + >>> plt.show() + + The numerical Laplace transform has an issue on the right, which is due to + the circular nature of the FFTLog integral transform. The effect is + mitigated by computing a biased transform with ``q = 0.5``. Good values of + the bias parameter ``q`` depend on the shape of the input function. + + >>> k, ak = fftl(u_laplace, r, ar, q=0.5) + >>> + >>> plt.loglog(k, ak) + >>> plt.loglog(k, lt, ':') + >>> plt.xlabel('$k$') + >>> plt.ylabel('$L[\\tanh](k)$') + >>> plt.show() + + ''' + + if np.ndim(r) != 1: + raise TypeError('r must be 1d array') + if np.shape(ar)[-1] != len(r): + raise TypeError('last axis of ar must agree with r') + + # inputs + n = len(r) + lnkr = np.log(kr) + + # make sure given r is logarithmic grid + if not np.allclose(r, np.geomspace(r[0], r[-1], n)): + raise ValueError('r it not a logarithmic grid') + + # log spacing + dlnr = (np.log(r[-1]) - np.log(r[0]))/(n-1) + + # low-ringing condition + if krgood: + y = PI/dlnr + um = np.exp(-1j*y*lnkr)*u(q + 1j*y, *args) + a = np.angle(um)/PI + lnkr = lnkr + dlnr*(a - np.round(a)) + + # transform factor + y = np.linspace(0, 2*PI*(n//2)/(n*dlnr), n//2+1) + um = np.exp(-1j*y*lnkr)*u(q + 1j*y, *args) + + # low-ringing kr should make last coefficient real + if krgood and not np.isclose(um[-1].imag, 0): + raise ValueError('unable to construct low-ringing transform, ' + 'try odd number of points or different q') + + # fix last coefficient to real when n is even + if not n & 1: + um.imag[-1] = 0 + + # bias input + if q != 0: + ar = ar*r**(-q) + + # set up k in log space + k = np.exp(lnkr)/r[::-1] + + # transform via real FFT + cm = np.fft.rfft(ar, axis=-1) + cm *= um + ak = np.fft.irfft(cm, n, axis=-1) + ak[..., :] = ak[..., ::-1] + + # debias output + ak /= k**(1+q) + + # output grid and transform + result = (k, ak) + + # derivative + if deriv: + cm *= -(1 + q + 1j*y) + dak = np.fft.irfft(cm, n, axis=-1) + dak[..., :] = dak[..., ::-1] + dak /= k**(1+q) + result = result + (dak,) + + # return chosen outputs + return result diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..71afcd3 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,3 @@ +[build-system] +requires = ['setuptools>=30.3.0', 'wheel'] +build-backend = 'setuptools.build_meta' diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..9218829 --- /dev/null +++ b/setup.cfg @@ -0,0 +1,33 @@ +[metadata] +name = fftl +version = attr:fftl.__version__ +maintainer = Nicolas Tessore +maintainer_email = n.tessore@ucl.ac.uk +description = generalised FFTLog +long_description = file: README.md +long_description_content_type = text/markdown +license = MIT +license_files = + LICENSE +url = https://github.com/ntessore/fftl +project_urls = + Documentation = https://fftl.readthedocs.io/ + Issues = https://github.com/ntessore/fftl/issues +classifiers = + Programming Language :: Python :: 3 + License :: OSI Approved :: MIT License + Operating System :: OS Independent + +[options] +python_requires = >=3.6 +install_requires = + numpy +py_modules = fftl + +[options.extras_require] +docs = + sphinx + furo + numpydoc + scipy + matplotlib diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..6068493 --- /dev/null +++ b/setup.py @@ -0,0 +1,3 @@ +from setuptools import setup + +setup()