From a3103ae7ba529b1a34c83ccbbdbad9abf53f264e Mon Sep 17 00:00:00 2001 From: Nicolas Tessore Date: Thu, 18 Mar 2021 13:05:40 +0000 Subject: [PATCH] initial commit --- LICENSE | 21 +++ MANIFEST.in | 3 + README.md | 20 +++ dctdlt.c | 107 ++++++++++++++ flt.pyx | 377 +++++++++++++++++++++++++++++++++++++++++++++++++ pyproject.toml | 2 + setup.cfg | 21 +++ setup.py | 6 + 8 files changed, 557 insertions(+) create mode 100644 LICENSE create mode 100644 MANIFEST.in create mode 100644 README.md create mode 100644 dctdlt.c create mode 100644 flt.pyx create mode 100644 pyproject.toml create mode 100644 setup.cfg create mode 100644 setup.py diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..59312e9 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2021 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/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..e7f8c58 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,3 @@ +include flt.pyx +include dctdlt.c +exclude flt.c diff --git a/README.md b/README.md new file mode 100644 index 0000000..fc4191b --- /dev/null +++ b/README.md @@ -0,0 +1,20 @@ + +flt +=== + +**fast Legendre transform** + +This is a minimal Python package for fast discrete Legendre transforms (DLTs). +The implementation uses a recursive version of the matrix relations by Alpert & +Rokhlin (1991) to compute the DLT via a discrete cosine transform (DCT). + +The package can be installed using pip: + + pip install flt + +For more information, please see the [documentation]. + +Current functionality covers the absolutely minimal use case. Please open an +issue on GitHub if you would like to see anything added. + +[documentation]: https://cltools.readthedocs.io/flt/ diff --git a/dctdlt.c b/dctdlt.c new file mode 100644 index 0000000..2ab2901 --- /dev/null +++ b/dctdlt.c @@ -0,0 +1,107 @@ +// dctdlt.c +// ======== +// discrete Legendre transform via DCT +// +// author: Nicolas Tessore +// license: MIT +// +// Synopsis +// -------- +// The `dctdlt` and `dltdct` functions convert the coefficients of a discrete +// cosine transform (DCT) to the coefficients of a discrete Legendre transform +// (DLT) and vice versa [1]. +// +// References +// ---------- +// [1] Alpert, B. K., & Rokhlin, V. (1991). A fast algorithm for the evaluation +// of Legendre expansions. SIAM Journal on Scientific and Statistical +// Computing, 12(1), 158-179. +// + +#define DCTDLT_VERSION 20210318L + + +// dctdlt +// ====== +// convert DCT coefficients to DLT coefficients +// +// Parameters +// ---------- +// n : unsigned int +// Length of the input array. +// dct : (n,) array of double +// Input DCT coefficients. +// dlt : (n,) array of double, output +// Output DLT coefficients. +// +void dctdlt(unsigned int n, const double* dct, double* dlt) +{ + double a, b; + unsigned int k, l; + + // first row + a = 1.; + b = a; + dlt[0] = 0.5*b*dct[0]; + for(k = 2; k < n; k += 2) + { + b *= (k-3.)/(k+1.); + dlt[0] += b*dct[k]; + } + + // remaining rows + for(l = 1; l < n; ++l) + { + a /= (1. - 0.5/l); + b = a; + dlt[l] = b*dct[l]; + for(k = l+2; k < n; k += 2) + { + b *= (k*(k+l-2.)*(k-l-3.))/((k-2.)*(k+l+1.)*(k-l)); + dlt[l] += b*dct[k]; + } + } +} + + +// dltdct +// ====== +// convert DLT coefficients to DCT coefficients +// +// Parameters +// ---------- +// n : unsigned int +// Length of the input array. +// dlt : (n,) array of double +// Input DLT coefficients. +// dct : (n,) array of double, output +// Output DCT coefficients. +// +void dltdct(unsigned int n, const double* dlt, double* dct) +{ + double a, b; + unsigned int k, l; + + // first row + a = 1.; + b = a; + dct[0] = b*dlt[0]; + for(l = 2; l < n; l += 2) + { + b *= ((l-1.)*(l-1.))/(l*l); + dct[0] += b*dlt[l]; + } + + // remaining rows + for(k = 1; k < n; ++k) + { + a *= (1. - 0.5/k); + b = a; + dct[k] = b*dlt[k]; + for(l = k+2; l < n; l += 2) + { + b *= ((l-k-1.)*(l+k-1.))/((l-k)*(l+k)); + dct[k] += b*dlt[l]; + } + } +} diff --git a/flt.pyx b/flt.pyx new file mode 100644 index 0000000..426455f --- /dev/null +++ b/flt.pyx @@ -0,0 +1,377 @@ +# flt: fast Legendre transform +# +# author: Nicolas Tessore +# license: MIT +# +# cython: language_level=3, boundscheck=False, embedsignature=True +# +''' + +Discrete Legendre Transform (:mod:`flt`) +======================================== + +This is a minimal Python package for fast discrete Legendre transforms (DLTs). +The implementation uses a recursive version of the matrix relations by Alpert & +Rokhlin (1991) to compute the DLT via a discrete cosine transform (DCT). + +The package can be installed using pip:: + + pip install flt + +Current functionality covers the absolutely minimal use case. Please open an +issue on GitHub if you would like to see anything added. + + +Reference/API +------------- + +.. autosummary:: + :toctree: api + :nosignatures: + + dlt + idlt + dltmtx + idltmtx + theta + +''' + +__all__ = [ + 'dlt', + 'idlt', + 'dltmtx', + 'idltmtx', + 'theta', +] + + +import numpy as np +from scipy.fft import dct, idct + +cdef extern from "dctdlt.c": + void dctdlt(unsigned int, const double*, double*) + void dltdct(unsigned int, const double*, double*) + + +def dlt(a, closed=False): + r'''discrete Legendre transform + + Takes a function in :math:`\mathtt{n}` points and returns its discrete + Legendre transform coefficients from :math:`0` to :math:`\mathtt{lmax} + = \mathtt{n}-1`. + + The function must be given at the points :math:`\cos(\theta)` returned by + :func:`flt.theta`. These can be distributed either over the open interval + :math:`\theta \in (0, \pi)`, or over the closed interval :math:`\theta \in + [0, \pi]`, in which case :math:`\theta_0 = 0`\ and :math:`\theta_{n-1} = + \pi`. + + Parameters + ---------- + a : (n,) array_like + Function values. + closed : bool, optional + Compute DLT over open (``closed=False``) or closed (``closed=True``) + interval. + + Returns + ------- + b : (n,) array_like + Legendre coefficients :math:`0` to :math:`\mathtt{lmax}`. + + See Also + -------- + flt.idlt : the inverse operation + flt.theta : compute the angles at which the function is evaluated + + Notes + ----- + The discrete Legendre transform takes a function :math:`f(\cos\theta)` over + the domain :math:`\theta \in [0, \pi]` and returns the coefficients of the + series + + .. math:: + + f(\cos\theta) = \sum_{l=0}^{l_{\max}} c_l \, P_l(\cos\theta) \;, + + where :math:`P_l(\cos\theta)` is a Legendre polynomial. + + The computation is done in two steps: + + First, the function is transformed to the coefficients of a discrete cosine + transform (DCT) using an inverse DCT-III for the open interval, or an + inverse DCT-I for the closed interval. + + Second, the DCT coefficients are transformed to the DLT coefficientsu using + a recursive version of the matrix relation given by [1]_. + + References + ---------- + .. [1] Alpert, B. K., & Rokhlin, V. (1991). A fast algorithm for the + evaluation of Legendre expansions. SIAM Journal on Scientific and + Statistical Computing, 12(1), 158-179. + + ''' + + # length n of the transform + if np.ndim(a) != 1: + raise TypeError('array must be 1d') + n = np.shape(a)[-1] + + # type of the DCT depends on open or closed interval + if closed: + dcttype = 1 + else: + dcttype = 3 + + # compute the DCT coefficients + a = idct(a, type=dcttype, axis=-1, norm=None) + + # this holds the DLT coefficients + b = np.empty(n, dtype=float) + + # these are memviews on a and b for C interop + cdef double[::1] a_ = a + cdef double[::1] b_ = b + + # transform DCT coefficients to DLT coefficients using C function + dctdlt(n, &a_[0], &b_[0]) + + # done + return b + + +def idlt(b, closed=False): + r'''inverse discrete Legendre transform + + Takes the :math:`\mathtt{n} = \mathtt{lmax}+1` coefficients of a DLT and + returns the corresponding function in :math:`\mathtt{n}` points. + + The function will be given at the points :math:`\cos(\theta)` returned by + :func:`flt.theta`. These can be distributed either over the open interval + :math:`\theta \in (0, \pi)`, or over the closed interval :math:`\theta \in + [0, \pi]`, in which case :math:`\theta_0 = 0` and :math:`\theta_{n-1} = + \pi`. + + Parameters + ---------- + b : (n,) array_like + DLT coefficients from :math:`0` to :math:`\mathtt{lmax}`. + closed : bool, optional + Compute function over open (``closed=False``) or closed + (``closed=True``) interval. + + Returns + ------- + a : (n,) array_like + Function values. + + See Also + -------- + flt.dlt : the forward operation + flt.theta : compute the angles at which the function is evaluated + + Notes + ----- + The inverse discrete Legendre transform returns a function + :math:`f(\cos\theta)` over the domain :math:`\theta \in [0, \pi]` given the + coefficients of the series + + .. math:: + + f(\cos\theta) = \sum_{l=0}^{l_{\max}} c_l \, P_l(\cos\theta) \;, + + where :math:`P_l(\cos\theta)` is a Legendre polynomial. + + The computation is done in two steps: + + First, the DLT coefficients are transformed to the coefficients of a + discrete cosine transform (DCT) using a recursive version of the matrix + relation given by [1]_. + + Second, the function values are computed using a DCT-III for the open + interval, or a DCT-I for the closed interval. + + References + ---------- + .. [1] Alpert, B. K., & Rokhlin, V. (1991). A fast algorithm for the + evaluation of Legendre expansions. SIAM Journal on Scientific and + Statistical Computing, 12(1), 158-179. + + ''' + + # length n of the transform + if np.ndim(b) != 1: + raise TypeError('array must be 1d') + n = np.shape(b)[-1] + + # type of the DCT depends on open or closed interval + if closed: + dcttype = 1 + else: + dcttype = 3 + + # this holds the DCT coefficients + a = np.empty(n, dtype=float) + + # these are memviews on a and b for C interop + cdef double[::1] a_ = a + cdef double[::1] b_ = b + + # transform DLT coefficients to DCT coefficients using C function + dltdct(n, &b_[0], &a_[0]) + + # perform the DCT + xi = dct(a, type=dcttype, axis=-1, norm=None) + + # done + return xi + + +def dltmtx(n, closed=False): + r'''discrete Legendre transform matrix + + Computes a matrix that performs the discrete Legendre transform + :func:`flt.dlt` when multiplied by a vector of function values. + + Parameters + ---------- + n : int + Length of the transform. + closed : bool, optional + Compute DLT over open (``closed=False``) or closed (``closed=True``) + interval. + + Returns + ------- + m : (n, n) array_like + Discrete Legendre transform matrix. + + See Also + -------- + flt.dlt : the equivalent operation + flt.theta : compute the angles at which the function is evaluated + + Notes + ----- + The discrete Legendre transform :func:`flt.dlt` performs the transformation + in place and does not compute the matrix :func:`flt.dltmtx`. + + ''' + + # type of the DCT depends on open or closed interval + if closed: + dcttype = 1 + else: + dcttype = 3 + + # compute the DCT matrix row by row (not column by column) + a = idct(np.eye(n), type=dcttype, axis=1, norm=None) + + # this holds the DLT matrix + b = np.empty((n, n), dtype=float) + + # these are memviews on a and b for C interop + cdef double[:, ::1] a_ = a + cdef double[:, ::1] b_ = b + + # transform DCT row to DLT row using C function + for i in range(n): + dctdlt(n, &a_[i, 0], &b_[i, 0]) + + # return transpose since we worked on rows, not columns + return b.T + + +def idltmtx(n, closed=False): + r'''inverse discrete Legendre transform matrix + + Computes a matrix that performs the inverse discrete Legendre transform + :func:`flt.idlt` when multiplied by a vector of coefficients. + + Parameters + ---------- + n : int + Length of the transform. + closed : bool, optional + Compute inverse DLT over open (``closed=False``) or closed + (``closed=True``) interval. + + Returns + ------- + m : (n, n) array_like + Inverse discrete Legendre transform matrix. + + See Also + -------- + flt.idlt : the equivalent operation + flt.theta : compute the angles at which the function is evaluated + + Notes + ----- + The inverse discrete Legendre transform :func:`flt.idlt` performs the + transformation in place and does not compute the matrix + :func:`flt.idltmtx`. + + ''' + + # type of the DCT depends on open or closed interval + if closed: + dcttype = 1 + else: + dcttype = 3 + + # this is the input matrix + b = np.eye(n, dtype=float) + + # this holds the DLT part of the matrix + a = np.empty((n, n), dtype=float) + + # these are memviews on a and b for C interop + cdef double[:, ::1] a_ = a + cdef double[:, ::1] b_ = b + + # transform DLT row to DCT row using C function + for i in range(n): + dltdct(n, &b_[i, 0], &a_[i, 0]) + + # multiply by DCT matrix + # return transpose since we worked on rows, not columns + return dct(a, type=dcttype, axis=1, norm=None).T + + +def theta(n, closed=False): + r'''compute angles for DLT function values + + Returns :math:`n` angles :math:`\theta_0, \ldots, \theta_{n-1}` at which + the function :math:`f(\cos\theta)` is evaluated in :func:`flt.dlt` or + :func:`flt.idlt`. + + The returned angles can be distributed either over the open interval + :math:`(0, \theta)`, or over the closed interval :math:`[0, \pi]`, in which + case :math:`\theta_0 = 0, \theta_{n-1} = \pi`. + + Parameters + ---------- + n : int + Number of nodes. + + Returns + ------- + theta : array_like (n,) + Angles in radians. + closed : bool, optional + Compute angles over open (``closed=False``) or closed (``closed=True``) + interval. + + ''' + + if closed: + t = np.linspace(0, np.pi, n, dtype=float) + else: + t = np.arange(n, dtype=float) + t += 0.5 + t *= np.pi/n + + return t diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..bc2be67 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,2 @@ +[build-system] +requires = ["setuptools", "wheel", "Cython"] diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..b99619f --- /dev/null +++ b/setup.cfg @@ -0,0 +1,21 @@ +[metadata] +name = flt +version = 2021.3.18 +maintainer = Nicolas Tessore +maintainer_email = n.tessore@ucl.ac.uk +description = fast Legendre transform +long_description = file: README.md +long_description_content_type = text/markdown +license = MIT +license_file = LICENSE +url = https://github.com/ntessore/flt +classifiers = + Programming Language :: Python :: 3 + License :: OSI Approved :: MIT License + Operating System :: OS Independent + +[options] +python_requires = >=3.6 +install_requires = + numpy + scipy diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..bc0a058 --- /dev/null +++ b/setup.py @@ -0,0 +1,6 @@ +from setuptools import setup +from Cython.Build import cythonize + +setup( + ext_modules=cythonize('flt.pyx'), +)